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A solution-adaptive,  algebraic  grid  generation  method  was  developed 
for  generating  grid  points  that  are  needed  by  finite -difference 
algorithms  to  obtain  numerical  solutions  to  complex  fluid  flow  problems. 
The  grid  generation  technique  developed  is  ideally  suited  for  use  in 
the  analysis  of  unsteady,  multidimensional  fluid  flow  problems  with 
deforming  spatial  domains  and  problems  with  disparate  length  scales  in 
different  parts  of  the  spatial  domain  which  can  move  about  in  an  un- 
predicable  manner. 

The  grid  generation  technique  developed  in  this  investigation  is 
based  on  variational  principles  and  does  not  require  solution  of  partial 
differential  equations  to  generate  grid  points.  The  method  will  auto- 
matically cluster  grid  points  in  regions  where  needed  to  resolve  all 
significant  length  scales  of  the  problem  being  analyzed.  In  addition, 


ix 


this  grid  generation  technique  will  ensure  that  the  grid  system  gener- 
ated is  smooth  and  nearly  orthogonal . 

The  orthogonality  of  the  grid  system  is  controlled  by  means  of 
parameters  introduced  to  ensure  that  the  deviation  from  orthogonality 
will  be  within  certain  limits,  fixed  a priori. 

To  demonstrate  the  usefulness  of  the  solution- adaptive , algebraic 
grid  generation  technique  developed  in  this  investigation,  the  method 
was  used  with  finite-difference  algorithms  to  analyze  (1)  the  linearized 
one -dimensional  inviscid  Burgers'  equation,  (2)  the  quasi-linear , one- 
dimensional inviscid  Burgers'  equation,  (3)  the  two-dimensional  un- 
steady problem  of  plane  wave  propagation,  and  (4)  the  propagation  of 
concentric,  circular  waves  in  the  radial  direction.  Also,  the  method 
was  used  in  conjunction  with  an  algebraic  grid  generation  method 
known  as  Two-Boundary  Technique  to  generate  a grid  network  in  a 
convergent -divergent -shaped  spatial  domain. 

Solutions  were  obtained  for  several  problems  and  compared  with  known 
solutions.  These  comparisons  indicated  that  the  solution- adaptive , al- 
gebraic grid  generation  method  developed  in  this  investigation  is  useful 
and  can,  efficiently,  generate  grid  points  required  by  finite -difference 
algorithms  to  obtain  numerical  solutions  for  fluid  flow  problems. 
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CHAPTER  I 


INTRODUCTION 
1 . 1 Background 

Computational  fluid  dynamics  has  experienced  tremendous  advances 
during  the  past  two  decades  and  has  emerged  as  one  of  the  most  powerful 
tools  for  analyzing  practical  fluid  flow  problems.  One  of  the  most 
important  factors  contributing  to  the  success  of  computational  fluid 
dynamics  as  a tool  for  analyzing  fluid  flow  problems  has  been  the 
advances  made  in  grid  generation  techniques.  The  evolution  of  grid 
generation  techniques  has  markedly  expanded  the  scope  of  fluid  flow 
problems  that  can  be  analyzed  by  numerical  methods.  Grid  generation 
techniques  supported  these  advances  because  they  enabled  computational 
fluid  dynamics  to  analyze,  efficiently  and  accurately,  problems  with 
very  complex  geometries,  problems  with  deforming  spatial  domain  as  well 
as  problems  with  highly  disparate  length  scales  in  different  parts  of 
the  domain.  Grid  generation  techniques  in  conjunction  with  finite - 
difference  methods  constitute  one  of  the  most  powerful  tools  for 
analyzing  practical  fluid  flow  problems  [1] . 

Before  the  advent  of  grid  generation  techniques,  complicated  in- 
terpolation schemes  were  needed  to  implement  boundary  conditions  when 
finite -difference  methods  were  used  to  analyze  fluid  flow  problems  with 
irregularly  shaped  spatial  domains.  At  that  time,  those  numerical 
methods  apparently  were  not  able  to  solve  fluid  flow  problems  with 
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complex- shaped  boundaries  or  problems  with  disparate  length  scales 
in  different  parts  of  spatial  domain  that  can  move  about.  The  develop- 
ment of  grid  generation  techniques  removed  those  limitations , gave  new 
impetus  to  numerical  methods  of  analysis,  and  made  finite -difference 
methods  versatile  and  accurate  in  analyzing  fluid  dynamics  problems  with 
irregularly  shaped  and/or  deforming  spatial  domains  [2], 

Grid  generation  is  concerned  with  the  mapping  of  optimally  dis- 
tributed grid  points  which  may  be  non-equally  spaced  and/or  moving  in 
the  physical  domain  onto  a computational  domain  where  all  of  the  grid 
points  are  uniformly  distributed  and  stationary.  This  invariably 
involves  the  generation  of  a coordinate  system  known  as  the  boundary- 
fitted  coordinate  system  for  the  computational  domain.  The  boundary- 
fitted  coordinate  system  is  said  to  be  generated  once  the  coordinate 
transformation  between  the  boundary- fitted  coordinate  system  and  the 
coordinate  system  of  the  physical  domain  is  known.  The  boundary- fitted 
coordinate  system  is  a curvilinear  coordinate  system  in  which  one  set  of 
coordinate  lines  or  surfaces,  analytically,  always  coincides  with  the 
boundaries  of  the  physical  domain  regardless  of  its  shape  or  motion  [3]. 

Once  the  boundary- fitted  coordinate  system  has  been  established,  the 
governing  equations  are  transformed  so  that  the  coordinates  of  the 
boundary- fitted  coordinate  system  become  the  independent  variables 
instead  of  the  coordinates  of  the  physical  domain.  Afterwards,  the 
finite-difference  equations  needed  to  obtain  numerical  solutions  are 
derived  from  these  transformed  governing  equations.  The  algebraic  equa- 
tions thus  derived  are  analyzed  in  the  computational  domain  where  all  of 
the  grid  points  are  stationary  and  uniformly  distributed.  Performing 
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calculations  in  this  computational  domain  is  important  because  it 
permits  the  development  of  general  numerical  codes  that  can  be  applied 
to  solve  many  different  fluid  flow  problems  with  different  geometries  by 
simply  changing  the  boundary- fitted  coordinate  system  and  the  boundary 
conditions . 

At  this  point,  it  is  important  to  note  that  boundary- fitted  coordi- 
nate systems  should  satisfy  several  properties  as  described  below. 

First,  the  coordinate  transformation  between  the  boundary- fitted 
coordinate  system  and  the  coordinate  system  of  the  physical  domain 
should  be  continuous  and  differentiable  with  a Jacobian  that  is  never 
zero  so  that  the  mapping  will  be  one  to  one. 

Second,  in  order  to  facilitate  implementation  of  boundary 
conditions,  one  set  of  grid  points  should  always  coincide  with  the 
boundaries  of  the  physical  domain  regardless  of  its  shape  or  motion. 
That  is  why  one  set  of  coordinate  lines  or  surface  of  the  boundary- 
fitted  coordinate  system  should  always  coincide  with  the  boundaries  of 
the  physical  domain. 

Third,  for  computational  efficiency,  the  number  of  grid  points  used 
should  be  kept  to  the  minimum  required  to  resolve  accurately  all 
significant  length  scales  of  the  problem.  Thus,  in  order  to  achieve  the 
same  accuracy,  grid  points  are  often  moving  and  nonuniformly  distributed 
within  the  physical  domain  with  more  grid  points  in  those  regions  where 
higher  resolution  is  needed.  The  process  of  concentrating  more  grid 
points  in  certain  regions  is  referred  to  as  the  stretching  of  the 
boundary- fitted  coordinate  system. 

Fourth,  when  grid  points  in  the  physical  domain  are  nonuniformly 
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distributed  because  of  stretching,  the  boundary- fitted  coordinate  system 
should  be  smooth;  i.e.,  the  spacing  between  grid  points  should  change 
gradually  from  a region  where  grid  spacings  are  small  to  a region  where 
grid  spacings  are  large.  Typically,  the  ratio  of  two  adjacent  grid 
spacings  should  be  approximately  between  0.5  and  2 [3]. 

Fifth,  the  boundary- fitted  coordinate  system  should  be  such  that  the 
grid  system  in  the  physical  domain  is  nearly  orthogonal.  This  means 
that,  the  angle  between  intersecting  grid  lines,  at  the  interior  grid 
points  should  be,  approximately,  between  45  and  135  degrees.  Mastin 
[4] , showed  that  the  factor  which  causes  an  increase  in  truncation  error 
due  to  departure  from  orthogonality  is  directly  proportional  to  the 
inverse  of  the  sine  of  the  angle  between  intersecting  grid  lines  in  the 
physical  domain.  Thus,  some  deviation  from  orthogonality  has  negli- 
gible effects  on  truncation  errors  and  many  examples  were  given  in  the 
reference  cited  to  demonstrate  this.  However,  grid  lines  should  always 
be,  approximately,  orthogonal  to  the  boundaries  of  the  physical  domain 
to  facilitate  implementation  of  derivative  boundary  conditions. 

The  ideal  boundary- fitted  coordinate  system  is  one  that  possesses 
all  of  the  five  properties  discussed  above.  In  addition,  the 
distribution  of  grid  points  in  the  physical  domain  should  be  coupled 
dynamically  to  the  numerical  solution  as  it  is  being  computed  so  that 
grid  points  in  the  physical  domain  will  automatically  cluster  in  regions 
where  they  are  needed  and  when  they  are  needed.  Such  a grid  system  is 
referred  to  as  a solution -adaptive  grid  system. 

For  many  important  fluid  flow  problems,  the  positions  of  those 
regions  requiring  a higher  concentration  of  grid  points  are  known  a 
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priori  and  a satisfactory  boundary- fitted  coordinate  system  can  be 
established  by  means  of  relatively  simple  techniques.  However,  for 
other  problems  in  which  those  regions  move  about  in  an  unpredictable 
manner,  a solution- adaptive  grid  generation  technique  is  needed. 
Solution- adaptive , grid  generation  techniques  generally  utilize 
variational  principles  to  satisfy  those  requirements  of  stretching, 
smoothness,  and  orthogonality  of  the  boundary- fitted  coordinate  system 
as  it  is  being  generated. 

All  solution-adaptive,  grid  generation  methods  developed  so  far  can 
be  classified  into  two  major  categories:  (1)  algebraic  grid  generation 
methods  and  (2)  differential  equation  grid  generation  methods.  Alge- 
braic grid  generation  methods  use  interpolation  techniques  to  generate 
boundary- fitted  coordinate  systems  and  can  provide  precise  controls  over 
the  distribution  of  grid  points  in  the  physical  domain  by  means  of 
stretching  functions.  Also,  they  are  highly  efficient  because  they  do 
not  need  to  solve  partial  differential  equations  to  generate  boundary- 
fitted  coordinate  systems.  Disadvantages  of  algebraic  grid  generation 
methods  include  difficulties  in  controlling  smoothness  and  ortho- 
gonality. 

Differential  equation  grid  generation  methods  generate  boundary- 
fitted  coordinate  systems  by  solving  systems  of  partial  differential 
equations.  Since  those  systems  of  partial  differential  equations  must 
be  solved  numerically,  the  differential  equation  grid  generation  methods 
are  less  efficient  than  algebraic  grid  generation  methods . For  three 
dimensional,  solution- adaptive  grid  systems  in  which  a different  grid 
system  must  be  generated  at  each  time  step,  differential  equation 
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grid  generation  methods  may  not  be  feasible,  because,  in  this  case  they 
become  computationally  very  expensive.  Most  solution- adaptive , grid 
generation  methods  developed  so  far  utilize  variational  principles  while 
constructing  boundary- fitted  coordinate  systems.  Variational  principles 
have  been  employed  to  optimize  those  requirements  of  stretching, 
smoothness  and  orthogonality  of  the  boundary- fitted  coordinate  system. 
Variational  techniques  utilized  while  generating  boundary-fitted 
coordinate  systems  can  be  used  with  either  differential  equation  grid 
generation  methods  or  algebraic  grid  generation  methods.  The  method 
developed  in  this  investigation  was  a solution- adaptive , algebraic  grid 
generation  method  and  variational  principles  were  utilized. 

1.2  Literature  Survey 

Efforts  have  been  made  during  this  past  decade  to  develop  grid 
generation  techniques . All  grid  generation  techniques  developed  so  far 
can  be  divided  into  two  major  categories:  (1)  non-adaptive , grid 
generation  methods  and  (2)  solution- adaptive  grid  generation  methods. 
Non-adaptive  grid  generation  methods  are  designed  for  steady  problems 
and  problems  in  which  regions  of  sharp  gradients  are  known  a priori  so 
that  grid  points  can  be  clustered  correctly  to  improve  the  accuracy  of 
the  numerical  solution  and  at  the  same  time  reduce  the  number  of  grid 
points  needed  to  solve  the  problem.  Solution- adaptive  grid  generation 
methods  are  designed  for  unsteady  fluid  flow  problems  with  deforming 
spatial  domain  and/or  problems  with  disparate  length  scales  in  different 
parts  of  the  spatial  domain  which  can  move  about  in  an  unpredictable 
manner.  As  noted  earlier,  solution-adaptive  grid  generation  methods  can 
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be  further  divided  into  the  following  categories:  (1)  algebraic  grid 
generation  methods  and  (2)  differential  equation  grid  generation 
methods . 

An  authoritative  and  comprehensive  review  of  grid  generation 
techniques  was  presented  by  Thompson,  Zahir,  Warsi  and  Mastin  [5]. 
Their  review  covered,  in  depth,  differential  equation  grid  generation 
methods,  algebraic  grid  generation  methods,  and  solution-adaptive  grid 
generation  techniques,  developed  before  1980.  Since  that  time  many 
important  investigations  have  been  done  in  grid  generation.  This 
literature  survey  discusses  the  advances  made  in  solution- adaptive  grid 
generation  methods  since  1980. 

Hindman  and  Spencer  [6]  developed  a solution- adaptive , grid 
generation  scheme  based  on  elliptic  partial  differential  equations  with 
a time -dependent  control  function.  The  time -dependent  control  function 
dynamically  couples  the  numerical  solution  and  the  speed  at  which  each 
grid  point  should  travel  to  solve  the  physics  of  the  problem  being 
investigated.  A new  boundary- fitted  coordinate  system  is  generated  at 
each  time  step  based  on  the  computed  grid  speeds.  This  method  suffers 
from  cumulative  errors  in  the  computed  grid  speeds  so  that  grid  point 
locations,  after  some  time  interval,  need  to  be  readjusted  by  using  a 
non-adaptive , grid  generation  technique.  Although  their  method  is 
general  and  can  be  applied  to  analyze  multidimensional  problems,  it  only 
has  been  applied  to  analyze  one -dimensional  problems.  This  is  because 
their  method  is  computationally  expensive  and  it  is  difficult  to  con- 
struct relationships  between  the  numerical  solution  and  the  grid  speed. 

Saltzman  and  Brackbill  [7-10]  used  variational  principles  to  develop 
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a solution- adaptive , differential  equation  grid  generation  method.  Their 
method  is  based  on  an  idea  proposed  by  Winslow  [11].  Basically,  a 
functional  is  constructed  by  summing  three  integrals  representing 
stretching,  smoothness  and  orthogonality  of  the  boundary- fitted 
coordinate  system  in  the  physical  domain.  The  minimization  of  that 
functional  results  in  a set  of  elliptic  partial  differential  equations 
whose  solution  gives  the  boundary- fitted  coordinate  system.  The 
elliptic  nature  of  the  partial  differential  equations  representing  the 
coordinate  transformation  ensures  a one-to-one  mapping  and  generates  a 
smooth  and  nearly  orthogonal  grid  system.  Saltzman  and  Brackbill 
applied  their  scheme  to  solve  many  different  problems  including  an 
inviscid  supersonic  flow.  In  this  problem,  gradients  were  used  to 
determine  the  stretching  of  the  boundary- fitted  coordinate  system  in 
order  to  concentrate  grid  points  in  regions  with  shock  waves.  The 
disadvantage  of  this  method  is  that  it  is  computationally  expensive, 
since  elliptic  partial  differential  equations  must  be  solved.  This  is 
especially  true  for  unsteady,  three-dimensional  problems,  in  which  a 
different  grid  system  has  to  be  generated  at  each  time  step. 

Bell  and  Shubin  [12]  developed  a solution- adaptive  grid  generation 
technique  similar  to  that  developed  by  Saltzman  and  Brackbill,  but 
they  only  applied  their  method  to  analyze  one -dimensional  problems. 
Within  the  functional  to  be  minimized  by  variational  techniques,  they 
introduced  a spatial  derivative  term  and  a time  derivative  term  to 
control  the  position  of  the  grid  points  in  the  physical  domain  at  each 
time  step.  Unfortunately,  it  is  very  difficult  to  construct  the  time 
derivative  term  in  the  functional  to  be  minimized.  Bell  and  Shubin 
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concluded  after  several  applications  of  their  method  that,  in  comparison 
with  non-adaptive  grid  generations  techniques,  their  method  expended  a 
lot  of  work  for  only  a modest  increase  in  accuracy. 

Rai  and  Anderson  [13]  developed  a very  interesting  solution- adaptive 
grid  generation  method.  In  their  method,  each  grid  point  has 
associated  with  it  an  attraction  or  repulsion  depending  upon  the 
difference  between  the  error  of  the  numerical  solution  at  each  grid 
point  and  the  average  error  over  all  grid  points.  The  spacings  between 
grid  points  are  decreased  or  increased  in  proportion  to  those 
differences  and  the  movements  of  the  grid  points  are  attenuated  by  an 
inverse  power  of  the  grid  spacing  in  the  physical  domain.  The  magnitude 
of  the  first  or  second  derivative  of  the  solution  was  used  as  a measure 
of  the  error  in  the  numerical  solution.  In  this  method,  smoothness  and 
orthogonality  of  the  grid  system  were  not  controlled.  Therefore,  in 
order  to  prevent  oscillations  and  distortions  in  the  numerical 
solutions,  this  method  imposes  limitations  on  the  movement  of  the  grid 
points  in  the  physical  domain  [14]. 

In  many  fluid  flow  problems,  disparate  length  scales  in  the  solution 
occur  only  in  one  coordinate  direction.  For  those  problems  the  grid 
spacing  along  each  coordinate  direction  can  be  controlled  indepen- 
dently. Dwyer  et  al . [15]  have  developed  a method  to  control  the 
distribution  of  grid  points  along  one  coordinate  direction  in  the 
physical  domain  by  using  a positive  weighting  function.  Such  positive 
weighting  function  was  a linear  combination  of  first  and  second 
derivatives  of  the  solution  with  respect  to  the  chosen  coordinate.  The 
method  has  been  applied  successfully  to  solve  many  problems  in  heat 
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transfer  and  in  two-dimensional  flame  propagation  problems.  The 
solution- adaptive , grid  system  generated  by  this  method  was  found  to 
improve  the  accuracy  of  numerical  solutions  considerably  without  the 
need  to  increase  the  number  of  grid  points  or  the  order  of  accuracy  of 
the  finite -differences  equations.  Dwyer  et  al.  also  noted  that  their 
solution- adaptive , grid  generation  method  can  remove  oscillations  often 
encountered  in  numerical  solutions. 

Gnoffo  [16,17]  developed  a solution- adaptive , grid  generation  method 
based  on  tension  spring  analogy  to  redistribute  grid  points  along  one 
coordinate  direction.  In  this  method,  the  weighting  function  used  was 
gradients  of  the  numerical  solution  which  are  treated  as  spring 
constants  in  a system  linking  the  grid  points  over  the  entire  physical 
domain.  To  increase  the  number  of  grid  points  in  certain  regions,  the 
spring  constants  were  made  larger  in  those  regions,  while  the  potential 
energy  of  the  spring  system  was  minimized.  This  method  was  applied  with 
success  to  compute  solutions  to  hypersonic  viscous  flow  problems.  These 
applications  showed  that  adaptation  along  only  one  coordinate  direction 
was  sufficient  to  move  the  grid  points  toward  regions  where  they  are 
needed  the  most  to  resolve  the  physics  of  the  problem  being  analyzed. 

Nakahashi  and  Deiwert  [18,19]  used  the  approach  of  Gnoffo  to  develop 
a similar  solution-adaptive,  grid  generation  method.  They  extended  the 
spring  analogy  to  solve  multidimensional  fluid  flow  problems  by  applying 
successive  adaptation  in  each  coordinate  direction.  The  smoothness  was 
controlled  by  introducing  constraints  to  ensure  grid  spacings  do  not 
change  too  fast  at  each  time  step  and  the  orthogonality  was  controlled 
by  introducing  torsional  springs  between  grid  points  in  the  physical 
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domain.  They  applied  their  method  to  generate  a grid  system  in  which 
grid  points  were  adapted  independently  in  each  coordinate  direction  to 
study  transonic  flows  past  an  airfoil.  The  density  gradient  was  used  to 
adapt  the  grid  system  in  the  streamwise  direction  for  the  purpose  of 
resolving  shock  waves.  In  the  direction  normal  to  the  airfoil  surface, 
the  component  of  momentum  in  that  direction  was  used  to  adapt  the  grid 
system  in  order  to  resolve  the  physics  in  the  boundary  layer  region.  The 
advantage  of  this  method  is  that  the  grid  system  is  adapted 
independently  in  each  coordinate  direction  and  this  permits  the  use  of 
marching  techniques.  Thus,  this  method  is  much  more  efficient  than 
methods  requiring  solution  of  elliptic  partial  differential  equations 
for  generating  grid  systems. 

Smooke  and  Koszykowski  [20]  developed  an  adaptive  grid  generation 
method  for  solving  one -dimensional  initial -boundary  value  problems  in 
which  the  number  of  grid  points  was  adjusted  to  equidistribute  a 
positive  weighting  function.  The  smoothness  was  achieved  by  imposing 
constraints  such  that  the  ratio  between  two  adjacent  grid  spacings  was 
less  than  a constant,  given  a priori,  and  greater  than  the  inverse  of 
that  constant.  The  method  required  interpolation  of  the  solution 
because  the  number  of  grid  points  are  not  the  same  at  each  time  level 
and,  sometimes,  this  was  found  to  cause  problems.  The  method  was 
applied  successfully  to  solve  several  combustion  problems  [21-22] . 

1.3  Objectives 

Even  though  considerable  progress  has  been  made  in  grid  generation, 
much  work  remains  to  be  done  in  developing  those  techniques.  This  is 
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especially  true  for  fluid  flow  problems  with  deforming  spatial  domain 
and  problems  with  highly  disparate  length  scales  in  different  part 
of  the  spatial  domain  that  can  move  about  in  a unpredictable  manner. 
The  major  objective  of  this  investigation  was  to  develop  a versatile  and 
efficient  solution- adaptive  grid  generation  technique  that  can  be  used 
in  the  analysis  of  such  fluid  flow  problems  using  finite-difference 
algorithms . 

The  grid  generation  method  developed  in  this  investigation  differs 
from  those  methods  developed  by  previous  investigators  in  the  following 
ways : 

1.  A mathematical  formulation  using  variational  principles  was 
applied  to  develop  a new  technique  for  generating  stretching 
functions.  The  stretching  functions  generated  by  using  this 
technique  couple  the  distribution  of  grid  points  in  the  physical 
domain  to  the  numerical  solution  as  it  is  being  computed. 

2.  A new  technique  was  developed  based  on  the  stretching  function 
mentioned  above  and  any  algebraic,  non-adaptive  grid  generation 
technique  (e.g.,  the  Two-Boundary  Technique)  to  create  a solution 
adaptive,  algebraic  grid  generation  method. 

3.  A new  technique  was  developed  to  ensure  that  the  grid  system  in 

the  physical  domain  is  smooth  and  that  the  grid  lines  do  not 
overlap.  Here,  two  new  concepts  were  introduced:  (1)  the 

weighting  function  was  controlled  by  a smoothness  factor  and  (2) 
a minimum  grid  spacing  specified  a priori  was  used  to  control  the 
nonoverlapping  of  the  grid  lines  in  those  problems  with  length 
scales  which  continuously  get  shorter  as  time  progresses. 
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4.  A new  technique  was  developed  to  ensure  that  the  grid  system  in 
the  physical  domain  is  nearly  orthogonal.  When  the  method 
developed  was  applied  in  conjunction  with  the  Two-Boundary 
Technique,  an  analytical  approach  was  developed  to  determine  the 
coefficients  for  controlling  the  orthogonality  at  the  boundaries. 

The  solution-adaptive,  algebraic  grid  generation  method  developed  in 
this  investigation  is  described  in  detail  in  Chapters  II  to  V and  is 
presented  in  the  following  order.  First,  the  technique  used  to 

transform  the  governing  equations  from  the  coordinates  of  the  physical 
domain  to  the  boundary- fitted  coordinates  of  the  computational  domain 
is  described.  Next,  the  procedure  by  which  variational  principles  are 
used  to  develop  stretching  functions  which  couple  the  distribution  of 
grid  points  to  the  numerical  solution  as  it  is  being  computed  is 
presented.  Then,  a description  is  given  on  how  stretching  functions 
based  on  variational  principles  are  used  to  construct  a solution- 
adaptive  grid  generation  method  for  one -dimensional  problems.  After- 
wards, how  the  solution- adaptive  grid  generation  technique  for  one- 
dimensional problems  can  be  applied  to  multidimensional  problems  is 
described.  Finally,  the  solution-adaptive  grid  generation  method 
developed  is  used  with  the  Two -Boundary  Technique  to  generate  grid 
points  inside  an  arbitrarily  shaped  spatial  domain. 

Applications  of  the  grid  generation  method  developed  in  this  inves- 
tigation are  described  in  Chapter  VI  to  VIII  and  the  conclusions  as  well 
as  recommended  future  works  are  given  in  Chapter  IX. 


CHAPTER  II 


COORDINATE  TRANSFORMATION  AND  METRIC  COEFFICIENTS 


In  this  chapter,  the  technique  used  to  transform  the  governing 
equations  from  the  coordinates  of  the  physical  domain  to  the  boundary- 
fitted  coordinates  of  the  computational  domain  is  described.  Also,  the 
technique  for  determining  the  metric  coefficients  for  a two-dimensional 
coordinate  transformation  is  presented. 

Grid  generation  invariably  involves  the  development  of  a boundary- 
fitted  coordinate  system  which  is  said  to  be  determined  once  the  co- 
ordinate transformation  between  the  coordinate  system  of  the  physical 
domain  and  the  boundary- fitted  coordinate  system  of  the  computational 
domain  is  known.  If  x-y-t  is  the  coordinate  system  of  the  physical 
domain  and  £-17-7-  is  the  boundary- fitted  coordinate  system  of  the 
computational  domain,  then  the  required  coordinate  transformation  has 
the  following  form: 

X - X ( £ , If  , T ) 

y - y (£.»?,  O (2.1) 

t - t ( r ) 

Each  grid  point  in  the  physical  domain  must  correspond  to  only  one 
grid  point  in  the  computational  domain  and  vice  versa.  This  means  that 
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the  coordinate  transformation  between  the  coordinate  system  of  the 
physical  domain  and  the  coordinate  system  of  computational  domain  must 
be  one  to  one.  The  necessary  and  sufficient  condition  for  a mapping  to 
be  one  to  one  is  that  the  Jacobian  of  the  coordinate  transformation  be 
nonzero  [23] . For  the  mapping  between  the  two  coordinate  systems 
expressed  by  Eq.  (2.1),  the  Jacobian  is  given  by 


3 ( x , y , t ) 

j = 

3 ( £ , v , t ) 


x£  xr?  xr 

y»?  yr 

0 0 tr 


- tr  ( x£  y^  - xr?  y^  ) 


(2.2) 


where  the  subscripts  denote  partial  derivatives. 

Since  the  Jacobian  of  the  coordinate  transformation  cannot  be  zero, 
the  inverse  of  the  coordinate  transformation  exists  and  has  the 
following  form: 


£ « £ ( x , y , t ) 

V - V ( X , y , t ) (2.3) 

r - t ( t ) 

Figure  2.1  shows  the  transformation  of  the  physical  domain  in  a 
x-y-t  coordinate  system  onto  the  computational  domain  in  a t 

boundary- fitted  coordinate  system.  The  grid  points  in  the  computational 
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Figure  2.1  Transformation  of  (a)  a physical  domain  to 
(b)  a computational  domain 
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domain  are  defined  by  the  intersection  of  equally  spaced  coordinate 
lines.  After  the  locations  of  grid  points  in  the  computational  domain 
are  known,  the  locations  of  grid  points  in  the  physical  domain  can  be 
determined  by  using  Eq.  (2.1). 

Once  the  coordinate  transformation  between  the  coordinate  systems  of 
the  physical  and  the  computational  domains  has  been  determined,  the 
governing  equations  for  the  problem  to  be  analyzed  must  be  transformed 
so  that  £,  r?  and  r become  the  independent  variables  instead  of  x,  y and 
t.  The  transformation  of  the  governing  equations  is  accomplished  in  two 
steps.  First,  the  relationship  between  partial  derivatives  with  respect 
to  x,  y and  t and  partial  derivatives  with  respect  to  f , ?/  and  r must  be 
established.  Second,  the  partial  derivatives  with  respect  to  x,  y and  t 
appearing  in  the  governing  equations  must  be  replaced  by  the  partial 
derivatives  with  respect  to  £ , r\  and  r. 

Any  partial  derivative  of  a dependent  variable  f with  respect  to  x, 
y and  t can  be  expressed  in  terms  of  £ , r]  and  r by  using  the  coordinate 
transformation  given  by  Eq.  (2.1)  and  the  chain  rule  as  follows: 


■ 

" 

af 

£x 

lx 

0 

af 

dx 

a* 

df 

riy 

0 

af 

dy 

dn 

df 

it 

r>y 

df 

at 

Tt 

dr 

_ 

. 

(2.4) 


In  the  equation  above 


the  elements  of  the  matrix  are  known  as 
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metric  coefficients  of  the  coordinate  transformation.  The  metric 
coefficients  depend  only  on  the  nature  of  the  coordinate  transformation 
given  by  Eq.  (2.1)  and  do  not  depend  on  the  physics  of  the  problem  being 
investigated.  Equation  (2.4)  can  be  written  in  a vectorial  form  as 

-4  -4 

f - [A]  g where  [A]  represent  the  matrix  of  the  metric  coefficients  of 
the  coordinate  transformation. 

Similarly,  by  using  the  inverse  of  the  coordinate  transformation 
given  by  Eq.  (2.3),  the  following  relation  can  be  obtained: 


“ 

■ 

af 

y* 

0 

af 

ax 

af 

0 

af 

dr) 

xn 

« 

— 

ay 

af 

Yr 

af 

dr 

xr 

at 

. 

(2.5) 


which  in  a vectorial  form  is  expressed  as  g = [B]  f. 

■4  *+  — > -4 

Using  Eqs . (2.4)  and  (2.5),  namely,  f - [A]  g and  g - [B]  f, 

-1 

gives  [A]  - [B]  . Since  the  elements  of  matrix  [A]  should  be  identical 
to  the  elements  of  the  inverse  of  the  matrix  [B]  , the  metric  coeffi- 
cients can  be  written  as  follows: 


Yr,  Yr 

-l 

xr) 

XT 

o 

rt 

1 

<-t  J 

0 

rt 

J 


(2.6,7) 
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(2.8,9) 


(2.10,11) 


x£  x»? 

ye  yv 


1 


t-T 


(2.12) 


Substituting  Eqs . (2.6)  through  (2.12)  into  Eq.  (2.4)  gives  the 
desired  relationships  between  the  partial  derivatives  with  respect  to  x, 
y and  t and  the  partial  derivatives  with  respect  to  £ , rj  and  r.  These 
relationships  are  used  to  replace  the  partial  derivatives  with  respect 
to  x,  y,  and  t appearing  in  the  governing  equations  by  partial  deriva- 
tives with  respect  to  £ , rj  and  r.  After  the  governing  equations  have 
been  expressed  in  terms  of  the  boundary- fitted  coordinates,  the 
governing  equations  are  transformed  to  finite-difference  equations  which 
are  used  to  analyze  the  problem  in  the  computational  domain. 

Here,  it  is  noted  that  the  metric  coefficients  given  by  Eqs.  (2.6) 
through  (2.12)  are  only  valid  for  unsteady,  two-dimensional  problems. 
However,  the  procedure  used  to  derive  those  equations  can  be  applied  in 
a straightforward  manner  to  derive  the  metric  coefficients  for  unsteady, 
three-dimensional  problems  [3]. 


CHAPTER  III 


THE  SAAGG  METHOD  FOR  ONE -DIMENSIONAL  PROBLEMS  / 

In  this  chapter,  a solution- adaptive , algebraic  grid  generation 
(SAAGG)  method  for  one -dimensional  problems  is  developed.  The  method  is 
presented  in  the  following  sequence.  First,  a new  class  of  stretching 
function  based  on  variational  principle  is  presented.  Second,  how  the 
stretching  function,  thus  developed,  is  used  to  create  a solution- adap- 
tive grid  generation  method  is  described.  Third,  a smoothness  factor 
is  defined  to  control  the  smoothness  of  the  grid  system  being  generated. 
Finally,  the  algorithm  for  this  grid  generation  method  is  summarized. 

3,1  Variational  Principles  and  Stretching  Functions 
Stretching  functions  have  been  widely  used  in  grid  generation 
techniques  to  redistribute  grid  points  in  the  physical  domain  to  improve 
accuracy  of  numerical  solution  by  concentrating  grid  points  where  most 
needed  [24],  In  order  to  construct  appropriate  stretching  functions, 
it  is  necessary  to  know  qualitatively  the  physics  of  the  problem  being 
analyzed.  For  unsteady  fluid  flow  problems  in  which  disparate  length 
scales  in  different  parts  of  the  physical  domain  can  move  about  in  an 
unpredictable  manner,  the  required  qualitative  information  is  unknown 
and  it  is  difficult  to  construct  appropriate  stretching  functions.  For 
such  problems,  the  stretching  function  should  be  coupled  to  the 
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numerical  solution  being  computed  and  it  also  should  be  able  to 
redistribute  any  existing  grid  system  to  minimize  numerical  errors.  A 
stretching  function  with  such  features,  derived  from  a variational 
technique,  is  developed  in  this  chapter  for  one -dimensional  problems. 

In  order  to  generate  a stretching  function  based  on  a variational 
technique,  some  functional  to  be  minimized  must  be  formed  first.  In 
this  investigation,  that  functional  was  taken  to  be  the  integral  of  the 
product  of  some  property  of  the  grid  (G) , related  to  the  grid  spacing  in 
the  physical  domain,  and  a positive  weighting  function  (W) , related  to 
the  error  in  the  numerical  solution.  The  minimization  of  this  func- 
tional guarantees  that  the  grid  points  in  the  physical  domain  will  be 
distributed  so  that  the  grid  spacing  will  be  small  where  the  weighting 
function  is  large,  and  vice  versa.  The  reason  W was  squared  in  the 
functional  below  is  to  ensure  that  the  weighting  function  is  always 
positive.  Thus,  the  functional  to  be  minimized  is  expressed  by 


I - f W (x)  G(x)  d£ 


(3.1.1) 


In  the  equation  above,  £ is  the  coordinate  of  computational  domain,  x 
is  the  coordinate  of  the  physical  domain  before  stretching,  and  x is  the 
coordinate  of  the  physical  domain  after  stretching.  Both  x and  x are 
functions  of  £.  The  property  of  grid  (G)  is  defined  below  and  the 
weighting  function  (W)  is  defined  in  Section  3.2.1. 

The  property  of  the  grid  (G)  in  the  functional  to  be  minimized  was 
chosen  to  be  the  square  of  the  grid  spacing  in  the  physical  domain  after 
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stretching.  If  the  grid  spacing  in  the  computational  domain  is  taken  to 
be  unity  (i.e.  A£-l) , then  the  grid  spacing  in  the  physical  domain  (Ax) 
can  be  expressed  as  follows : 

Ax 

= or  Ax  « Xtf  (3.1.2) 

A? 


Thus,  the  property  of  the  grid  (G)  can  be  written  as  follows: 


(3.1.3) 


The  stretching  function  (K)  to  be  developed  in  this  investigation 
relates  the  coordinates  of  the  physical  domain  before  and  after 
stretching  as  follows: 


x - K(x)  x (3.1.4) 

By  using  Eqs.  (3.1.3)  and  (3.1.4),  the  functional  to  be  minimized  by 
the  variational  technique  expressed  by  Eq.  (3.1.1)  can  now  be  written  as 
shown  below: 

. 2 
I - f (W(x)  [K(x)x] £ ) d£ 


(3.1.5) 
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The  partial  derivative  term  appearing  in  Eq.  (3.1.5)  is  obtained  by 
differentiation  as  follows: 


[K(x)x] £ 


3x  3K 

K(x)  + x 


3x 

K(x)  

a? 


dK  3x 

+ x 

dx  3£ 


dK 

[ K(x)  + x ] x^ 

dx 


Substitution  of  Eq.  (3.1.6)  into  Eq.  (3.1.5)  gives 


. 2 
I - f tW(x)  M(x)  x^  ] d£ 


where 


dK 

M(x)  - K(x)  + x 

dx 


(3.1.6) 


(3.1.7) 


(3.1.8) 


Equation  (3.1.7)  can  be  written  more  compactly  by  letting 


N(x)  - W(x)  M(x) 


(3.1.9) 
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so  that  the  functional  to  be  minimized  becomes 


I “ S t N(x)  X£  ] d£ 


(3.1.10) 


The  problem  of  minimizing  the  integral  above  can  be  solved  by  a 
variational  technique  and  the  solution  is  the  Euler -Lagrange  equation 
associated  with  that  integral.  The  Euler -Lagrange  equation,  derived  in 
Appendix  A,  is  as  follows: 


where 


and 


3F 


3x 


3F 

( ) 

3x£ 


0 


2 

F - [ N(x)  x^] 


3F  2 dN 

= 2 N (x)  x^  

3x  dx 


(3.1.11) 


(3.1.12) 


(3.1.13) 


3F 

Sx^ 


2 

2 N(x)  X| 


(3.1.14) 


Substituting  Eqs . (3.1.13)  and  (3.1.14)  into  Eq.  (3.1.11)  gives 


2 dN  d 2 

N(x)  x^  - [ N(x)  xc  ] 

dx  d^ 


0 


(3.1.15) 
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which  after  simplification  yields  the  following  ordinary  differential 
equation: 


d 

[ N(x)  X£  ] - 0 (3.1.16) 

dx 


By  integrating  the  above  equation  and  using  Eqs.  (3.1.8)  and 
(3.1.9),  Eq.  (3.1.16)  becomes 


^ £x 

[ K x ] - C2  (3.1.17) 

dx  W(x) 


where  C2  is  a constant  of  integration.  Equation  (3.1.17)  can  be 
integrated  again  to  yield 


1 

K(x) [ Ci  + C2  J dx  ] (3.1.18) 

x W(x) 

producing  another  constant  of  integration,  C]_. 

The  constants  of  integration,  and  C2,  can  be  evaluated  by 

requiring  that  the  coordinates  before  and  after  stretching  have  the  same 
value  at  each  boundary  of  the  physical  domain.  By  assuming  the 
boundaries  of  the  physical  domain  to  be  at  and  L2 , Eq.  (3.1.4) 
implies  that 
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K(x=L!)  - 1 


(3.1.19) 


and 


K(x=L2)  - 1 


(3.1.20) 


By  using  Eqs.  (3.1.19)  and  (3.1.20),  Eq.  (3.1.18)  yields  the 
following  system  of  equations: 


L1  “ Ci  + C2  [ J 


L2  = C]_  + C2  [ / 


W(x) 


dx  ] 


x=Li 


dx  ] 

W(x)  x-l2 


(3.1.21) 


(3.1.22) 


Solving  the  above  system  of  equations  for  the  constants,  and  C2 , 

yields 


and 


[ J 


Cl  - Li  - (L2-LX) 


W(x) 


dx  ] 


x-Li 


;Lz— 

Ll 


dx 


W(x) 


(3.1.23) 


C2  - 


l2  - L1 
L2  Cx 

f dx 

Li  W(x) 


(3.1.24) 
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Finally,  substituting  Eqs . (3.1.23)  and  (3.1.24)  into  Eq.  (3.1.18) 
gives  the  desired  expression  for  the  stretching  function  which  is 


dx 

.Lx  W(x) 

L1  + (l2  - Ll>  

L2  Cx 

dx 

. L]_  W(x) 

K(x)  - (3.1.25) 

x 


The  integrals  in  Eq.  (3.1.25)  can  be  evaluated  numerically  or 
analytically  to  determine  the  value  of  K(x)  at  each  grid  point.  Once 
K(x)  has  been  obtained,  it  is  substituted  into  Eq.  (3.1.4)  to  give  the 
location  of  grid  points  in  the  physical  domain  after  stretching. 

Application  of  this  technique  for  generating  stretching  functions  is 
presented  in  Chapter  VI  where  several  examples  are  given. 

3,2  A Solution-Adaptive  Grid  System 
The  stretching  function  presented  in  the  previous  section  can  couple 
the  distribution  of  grid  points  to  the  numerical  solution  of  the  problem 
that  is  being  analyzed  through  the  weighting  function.  Thus,  these 
stretching  functions  can  be  employed  to  develop  a solution-adaptive, 
grid  generation  method  to  be  used  for  those  problems  with  disparate 
length  scales  in  different  parts  of  the  spatial  domain  which  can  move 
about  in  an  unpredictable  manner.  This  solution- adaptive , grid 
generation  method  has  the  advantage  of  being  an  algebraic  method  because 
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partial  differential  equations  do  not  need  to  be  solved  in  order  to 
determine  the  locations  of  the  grid  points . 

As  noted  earlier,  in  order  for  a grid  generation  method  to  be 
useful,  it  should  satisfy  constraints  on  smoothness  and  orthogonality  in 
addition  to  stretching.  For  one -dimensional  problems,  only  stretching 
and  smoothness  need  to  be  considered.  In  this  section,  stretching  and 
smoothness  of  a one -dimensional  grid  system  are  discussed.  Stretching, 
smoothness  and  orthogonality  for  multidimensional  grid  systems  are 
discussed  in  the  next  chapter. 

The  stretching  function  represented  by  Eq.  (3.1.25)  can  be  substitu- 
ted into  Eq.  (2.1)  to  give  the  coordinate  transformation  necessary  to 
generate  grid  points  in  the  physical  domain.  By  noting  that 


1 

dx  - xf  d( d£  (3.2.1) 

^x 


the  relationship  for  the  coordinate  transformation  can  be  represented  by 
the  following  equation: 


x 


L1  + (l2  * Ll) 


f*  1 

d£ 

Jl  W(x) 


fIL  1 

d? 

. 1 W(x) 


(3.2.2) 
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where 


n+1 

x - x(  $ . r ) (3.2.3) 

and 

n+1 

x - x(  { , r ) (3.2.4) 


In  the  above  equations,  x is  the  coordinate  of  the  grid  points  in  the 
physical  domain  at  time  level  n+1  before  stretching;  x is  the  coordinate 
of  grid  points  at  time  level  n+1  after  stretching;  Lq  and  L2  are  the 
boundaries  of  the  spatial  domain;  and  1 and  IL  are  the  boundaries  of  the 
computational  domain. 

The  locations  of  the  grid  points  in  the  physical  domain,  determined 
by  Eq.  (3.2.2),  can  be  coupled  to  the  numerical  solution  as  it  is  being 
computed  through  the  weighting  function  (W)  . Such  weighting  function 
which  makes  the  method  to  be  a solution- adaptive  grid  generation  method 
is  presented  in  the  next  section. 


3,2.1  Weighting  Function 

The  weighting  function  appearing  in  Eq.  (3.2.2)  should  be  related  to 
the  numerical  solution  to  make  the  method  solution- adaptive . In  this 
investigation,  the  weighting  function  was  chosen  as  a function  of  the 
gradient  of  some  dependent  variable  of  the  problem  being  solved.  This 
is  because,  in  general,  numerical  errors  are  larger  in  those  regions  of 
the  physical  domain  where  the  gradients  are  larger  [25],  However,  let- 
ting the  weighting  function  equal  to  that  gradient  causes  two  problems. 
First,  when  the  gradient  is  zero,  Eq.  (3.2.2)  is  singular.  Second,  when 
the  gradient  is  very  large,  almost  all  of  the  grid  points  will  be 


30 


clustered  in  the  vicinity  of  the  large  gradient  so  that  the  grid  system 
will  not  be  smooth.  To  overcome  the  first  problem,  a positive  constant, 
in  this  investigation  chosen  to  be  unity,  was  added  to  the  absolute 
value  of  the  gradient  of  the  solution.  To  overcome  the  second  problem, 
the  constant  and  the  absolute  value  of  the  gradient  were  raised  to  some 
power  SF.  Thus,  the  weighting  function  used  in  this  investigation  has 
the  following  form: 


The  power  SF  is  referred  to  as  the  smoothness  factor  and  controls 
the  value  of  the  weighting  function  in  order  to  ensure  that  the  grid 
system  generated  is  smooth.  The  procedure  for  determining  the  smooth- 
ness factor  is  described  in  the  next  section. 

3,2.2  Smoothness  Factor 

The  spacing  between  two  adjacent  grid  points  in  the  physical  domain 
can  be  obtained  from  Eq.  (3.2.2)  as  follows: 


SF 


SF 


( 1 + I Ux  | ) 


(3.2.5) 


1 


Axf  = xi+1  - Xi  - (L2  - Lj_) 


TL 


(3.2.6) 


1 


d? 


SF 


1 


[ w ] 
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The  integrals  in  Eq.  (3.2.6)  can  be  evaluated  numerically  by  using 
the  backward- step  method  as  shown  below: 


Ax£  - (L2  - Li) 


SF 

1 


u>i 


SF 

IL  1 

2 

j-2  Wj 


(3.2.7) 


Equation  (3.2.7)  can  be  written  more  compactly  as  follows: 


Ax^  - 


L2  - LX 


SF 

IL 

2 

j-2 


(3.2.8) 


Here,  it  is  noted  that  higher-order  accurate  integration  methods 
could  have  been  used  instead  of  the  backward- step  method.  However,  for 
the  purpose  of  grid  generation,  the  backward- step  method  was  found  to 
give  results  as  good  as  those  obtained  by  using  Simpson's  rule. 

The  grid  spacings  in  the  physical  domain  can  be  no  greater  than 
^xmax>  and  no  smaller  than  Axm^n.  Both  Axmax  and  Axm^n  are  specified  a 
priori.  Hence,  since  is  a minimum  when  Ax^  is  a maximum  and  vice 
versa,  Eq.  (3.2.8)  gives  the  following  relations: 
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and 


IL 

2 


j-2 


SF 

“max  l2  ‘ L1 
“j  ^xmin 


IL 

S 

j-2 


SF 

“min  l2  * L1 
“j  Axmax 


(3.2.9) 


(3.2.10) 


Equations  (3.2.9)  and  (3.2.10)  can  be  rewritten  as  shown  below: 


SF  IL 
(wmax)  2 

j-2 


SF  IL 
<wmin)  2 

j-2 


0)4 


0)1 


SF 


l2  - L1 


Ax. 


mm 


SF 


L2  - LX 


Ax, 


max 


(3.2.11) 


(3.2.12) 


The  solution  of  the  above  system  of  equations  gives  the  following 
expression  for  the  smoothness  factor  (SF): 


Axmax 
In  

Axmin 

SF (3.2.13) 

“max 

In  


“min 
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The  above  expression  relates  SF  to  two  ratios,  namely  wmax/wm£n  and 
^xmax/^xmin-  A numerical  value  for  SF  can  be  obtained  once  these  two 
ratios  are  specified. 

The  ratio  WmaxA’min  is  determined  by  using  Eq . (3.2.5)  and  by 
examining  the  numerical  solution  (U)  being  computed  at  time  level  n,  as 
shown  below. 


wmax  ( 1 + max|Ux|) 

wmin  ( 1 + min | Ux | ) 


(3.2.14) 


In  the  above  equation,  max|ux|n  and  min|ux|n  denote  the  maximum  and  the 
minimum  absolute  values  of  Ux  in  the  physical  domain  at  time  level  n, 
respectively. 

Specification  of  the  ratio  Axmax/Axm^n  depends  on  the  problem  being 
analyzed.  If  the  maximum  gradient  of  the  problem  being  solved  does  not 
change  in  magnitude  as  time  progresses,  then  the  ratio  Axmax/Axmin  can 
be  set  equal  to  a constant  which  is  the  same  for  every  time  step;  i.e., 


Axmax 

= constant  (3.2.15) 

Axmin 

For  problems  in  which  the  maximum  gradient  of  the  problem  being 
solved  increases  as  time  progresses,  the  ratio  Axmax/Axmin  should  be 
changed  from  time  step  to  time  step.  Here,  that  ratio  is  specified  as 
follows : 
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n+1 


“max 

n+1 

n 

Axmax 

^xmax 

“min 

^xmin 

Axmin 

“max 

“min 

(3.2.16) 


where  (wmax)n+''"  and  (wniin)n+F  are  obtained  from  the  numerical  solution 
Un+1  by  using  a grid  system  obtained  at  time  level  n+1  before 
stretching. 

With  the  ratios  «niax/“min  and  ^max/Amin  specified  in  the  manner 
described  above,  Eq.  (3.2.13)  readily  gives  a numerical  value  for  SF. 
Once  SF  has  been  thus  determined,  the  location  of  each  grid  point  can  be 
determined  by  using  Eq.  (3.2.2)  as  follows: 


n+1 

xi  “ L1  + (^2  * Lx) 


[ *SF] 


n 


■IL 

Jl 


[ “SF1 


n 


(3.2.17) 


where  is  obtained  from  Eq.  (3.2.5). 

The  grid  generation  method  given  by  Eq.  (3.2.17)  may  produce  grid 
spacings  that  are  smaller  than  that  needed  to  obtain  solutions  of  the 
desired  accuracy.  In  order  to  prevent  this  situation,  the  smallest  grid 
spacing  allowed  (Axcrt)  must  be  specified  a priori.  The  precise  value 
of  Axcrt  depends  on  the  physics  of  the  problem  being  analyzed  and  the 
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accuracy  desired  in  the  computed  numerical  solution.  If  Axm^n  computed 
by  using  Eq.  (3.2.7)  is  less  than  Axcrt,  then  the  smoothness  factor  (SF) 
is  obtained  iteratively  from  the  following  nonlinear  equation: 


IL 

Axcrt  2 
j-2 


L2  - L]_ 

I [wmax]k+1| 


SF 


(3.2.18) 


where  k denotes  iteration.  The  above  equation  was  obtained  from  Eq. 
(3.2.9)  by  setting  Axmj_n  to  Axcrt. 

This  concludes  the  description  of  the  solution-adaptive  grid 
generation  method  for  one -dimensional  problems.  The  details  of  the 
implementation  of  this  method  are  described  in  the  next  section. 

3,3  Algorithm  for  One-Dimensional  Problems 

This  section  presents  the  algorithm  for  generating  solution-adaptive 
grid  systems  for  one -dimensional  problems  in  which  the  two  boundaries  of 
the  physical  domain  can  move  or  be  stationary.  The  algorithm  which  will 
map  the  grid  points  between  the  physical  domain  in  x-t  coordinate  system 
and  the  computational  domain  in  the  £— r coordinate  system  is  as  follows: 

(1)  Transform  the  governing  equations  from  the  coordinate  system  of 
the  physical  domain  to  the  coordinate  system  of  the  computa- 
tional domain. 

(2)  Apply  any  suitable  finite-difference  method  to  express  the 
governing  equations  as  a set  of  finite  difference  equations 
(FDEs) . 
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(3)  Specify  the  duration  of  interest  (i.e.  0<t<T)  and  set  the  time 
level  n to  zero  which  corresponds  to  time  t=0. 

(4)  Specify  the  number  of  grid  points  (IL)  to  be  employed.  IL 
depends  on  the  accuracy  desired  in  the  numerical  solutions  and 
on  the  computer  resources  available.  With  IL  specified,  Eq. 
(3.2.2)  implies  that  x-L^  corresponds  to  £-=l  and  x=L2  corre- 
sponds to  £=IL  as  expressed  by  the  following  equations: 


x (|-1)  - Li 


(3.3.1) 


x (£=IL)=  L2 


(3.3.2) 


(5)  Specify  the  locations  of  the  grid  points  in  the  computational 

domain.  As  noted  before,  the  grid  points  in  the  computational 

domain  are  always  equally  spaced  and  the  grid  spacing  is  unity 
(A£=l)  so  that  the  locations  of  the  grid  point  are  given  by 

£i  - i . i “ 1,2 IL  (3.3.3) 

(6)  Specify  the  locations  of  the  grid  points  in  the  physical  domain 

at  the  initial  time  level  (n=0) . The  initial  grid  system 

specified  may  consist  of  equally  spaced  or  nonequally  spaced 
grid  points.  If  the  grid  points  are  equally  spaced,  then  the 
locations  of  the  grid  point  are  given  by 
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L2  - LX 

*i  “ Lx  + (Cj.  - 1)  (3.3.4) 

IL  - 1 

i - 1,2 IL 

where  £x  is  given  by  Eq.  (3.3.3).  If  it  is  desired  to  use 
nonequally  spaced  grid  points  for  the  initial  grid  system 
because  of  the  nature  of  the  initial  conditions,  then  Eq. 
(3.2.17)  can  be  used  to  generate  such  a grid  system  as  shown 
below: 


xi  - Lx 


(3.3.5) 


n n n n 
xi  - Lx  + (L2  - Lx) 


1 1 

2 

n n n »-2  w(xm)SF 

« Lx  + (L2  - LX)  i-2,4 IL  (3.3.6) 

IL  1 

2 

m— 2 w(xm)SF 


where  is  given  by  Eq.  (3.3.3)  and 
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w(xm) 


1 + 


au(xra> 

dx 


* 1 + 


xm+l  * xm-l 


(3.3.7) 


In  the  above  equations,  U is  the  initial  condition;  xm  is  given 
by  Eq.  (3.3.4);  and  SF  is  determined  by  using  either  Eq. (3.2.13) 
or  Eq.  (3.2.18). 

(7)  Evaluate  the  metric  coefficients  numerically  from  the  following 
equations : 


2A£i 


n n 

xi+l  - xi-l 


and 


n+1  n 
2A£i  x±  - X£ 


n n 

xi+l  ‘ xi-l  Ati 


where 


(3.3.8) 


(3.3.9) 


n+1  n+1 

n+1  L2  " L1  n 

Xi  “ (x 

n n 

L2  - Li 


n n+1 

- Li)  + LX  (3.3.10) 


n 

In  the  above  equations  x^  is  obtained  in  step  6. 
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(8)  Use  the  FDEs  derived  in  step  2 to  obtain  a "tentative"  solution 
of  the  problem  at  the  next  level  by  using  the  grid  system 
generated  in  step  6. 

(9)  Determine  the  locations  of  the  grid  points  at  next  time  level 
by  using  Eqs . (3.3.6)  and  (3.3.7)  in  which  U(x)  is  the  solution 
obtained  in  step  8. 

n n+1 

(10)  Recalculate  the  metric  by  using  Eq.  (3.3.9)  with  x^  obtained 
in  step  9. 

(11)  Use  the  FDEs  with  the  metric  coefficients  given  by  Eqs.  (3.3.8) 
and  (3.3.9)  to  obtain  the  "corrected"  solution  at  the  next  time 
level . 

(12)  Replace  n by  n+1  and  calculate  the  time  at  that  time  level  by 
using  the  expression 

n+1  n+1 

t - S Atffi  (3.3.11) 

m-0 

n+1  n+1 

(13)  Repeat  steps  7 to  11  if  t < T.  If  t > T,  then  stop  the 
algorithm. 

Several  examples  illustrating  how  the  solution-adaptive  grid  genera- 
tion method  developed  in  this  chapter  can  be  applied  to  analyze  fluid 
flow  problems  are  presented  in  Chapter  VII. 


CHAPTER  IV 


THE  SAAGG  METHOD  FOR  2-D,  RECTANGULARLY  SHAPED  DOMAINS 

In  this  chapter,  a solution-adaptive,  algebraic  grid  generation 
(SAAGG)  method  for  two-dimensional,  rectangularly  shaped  spatial  domains 
is  presented.  This  method  is  based  on  the  one  - dimens ional , solution- 
adaptive  grid  generation  method  described  in  Chapter  III.  The  two- 
dimensional,  solution-adaptive  grid  generation  method  is  presented  in 
the  following  sequence.  First,  an  overview  of  the  method  is  de- 
scribed. Second,  the  techniques  used  to  control  the  smoothness  of  the 
grid  system  and  the  nonoverlapping  of  the  grid  lines  are  presented. 
Finally,  the  automatic  adjustments  for  controlling  the  orthogonality  of 
the  grid  system  are  defined. 

4.1  A Two-Dimensional.  Solution-Adaptive  Grid  System 
In  order  to  explain  the  solution-adaptive,  grid  generation  method 
for  two-dimensional,  rectangularly  shaped  spatial  domains,  consider 
mapping  the  physical  domain  shown  in  Fig.  4.1a  in  the  x-y-t  coordinate 
system  to  the  computational  domain  in  the  coordinate  system  shown 

in  Fig.  4.1b.  The  spatial  domain  shown  in  Fig  4.1a  can  undergo  recti- 
linear deformations.  Here,  the  coordinate  transformation  sought  between 
the  x-y-t  and  the  coordinate  systems  is  given  by  Eq.  (2.1)  and  is 

repeated  below  for  convenience. 
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(b) 

Figure  4.1  Correspondence  between  the  boundaries  of  the 
physical  domain  and  the  computational  domain. 

(a)  Rectangularly  shaped  physical  domain  (x-y-t) . 

(b)  Computational  domain  (£-r7~r)  . 
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x - x(£  , *? , t) 

y - yU.V.r)  (2.1) 

t - r 

Implementation  of  the  solution- adaptive , grid  generation  method  for 
two-dimensional,  rectangularly  shaped  spatial  domains  involves  the 
following  steps: 

(1)  Transform  the  governing  equations  from  the  coordinate  system  of 
the  physical  domain  to  the  coordinate  system  of  the  computa- 
tional domain. 

(2)  Apply  any  suitable  finite- difference  method  to  express  the 
transformed  governing  equations  as  a set  of  finite -difference 
equations  (FDEs) . 

(3)  Decide  the  correspondence  between  the  boundaries  of  the  physical 
domain  and  the  boundaries  of  the  computational  domain.  In  this 
investigation,  this  correspondence  is  shown  in  Fig.  4.1. 

(4)  Specify  the  number  of  grid  lines  IL  in  the  ^-direction  and  the 
number  of  grid  line  JL  in  the  ij- direction. 

(5)  Set  the  initial  time  level  n equal  to  zero  which  corresponds  to 
time  t*=0  and  specify  the  duration  of  interest  T (i.e.  0<t<T) . 

(6)  Specify  the  locations  of  the  grid  points  in  the  computational 
domain.  Here,  A£“l  and  Arj=l  so  that 

£i  “ i i-1,2 IL  (4.1.1) 

and 


r?j  - j 


j-1,2 JL 


(4.1.2) 
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(7)  Specify  the  locations  of  the  grid  points  in  the  physical  domain 
at  the  initial  time  level  (n-0)  . The  initial  grid  system  may 
have  grid  lines  that  are  equally  spaced  or  nonequally  spaced. 
If  it  is  desired  to  have  equally  spaced  grid  lines,  then  the  x- 
and  y-coordinates  of  the  grid  points  in  the  physical  domain  are 
given  by  the  following  expressions: 

n n 

n n ^2  - Li 

*i  - ^ + (£i  - 1)  (4.1.3) 

IL  - 1 

i - 1,2 IL 


,11 

n ,n  ^2  ' Lq 

yj  - Li  + (rn  - 1)  (4.1.4) 

JL  - 1 

j - 1.2 JL 

If  it  is  desired  to  have  nonequally  spaced  grid  lines,  because 
of  steep  gradients  appearing  in  the  initial  condition,  then  the 
x-coordinates  along  the  j ^ grid  line  are  given  by  Eqs . (3.3.5) 
to  (3.3.7).  In  those  equations  xm  is  given  by  Eq.  (4.1.3);  the 
smoothness  factor  is  taken  to  be  unity  (i.e.  SF=1) ; and  U(xm,xq) 
is  the  initial  condition  of  the  problem.  The  y-coordinates  of 
the  grid  points  are  obtained  by  using  the  procedure  just  de- 
scribed except  replace  x by  y,  i by  j,  Lq  by  l{ , L2  by  L2 , 
and  U(xm,yj)  by  U(xi,ym)  with  ym  given  by  Eq.  (4.1.4). 

(8)  Evaluate  the  metric  coefficients  of  the  coordinate  transforma- 
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tion  by  using  FDEs  obtained  from  Eqs . (2.6)  to  (2.12).  These 

FDEs  are  derived  by  using  second-order-accurate,  central- 

difference  formulas,  similar  to  that  shown  for  Eqs.  (3.3.8)  and 

n n 

(3.3.9).  The  coordinates  xjj  and  yjj  appearing  in  those  FDEs 

describe  the  locations  of  the  grid  points  at  time  level  n.  The 
n+1 

coordinate  xjj  appearing  in  those  FDEs  are  given  by  the  follow- 
ing equations : 


x 


n+1  n+1 

ij  “ Lij 


i-1 


(4.1.5) 


n+1 

{ij 


ALi 


n+1 


n n+1 

ASij  + Lij 


(4.1.6) 


n 


ALi 


i-2,3. . . ,IL 


where 


n IL  n 
ALj  - Z [ (x, 
m—z 


m.j 


n 2 n 

xm-l,i)  + <y, 


m,J 


n 2 1/2 

ym-l,j)  ] 


(4.1.7) 


AS  a 


n 


ij  ” t (xm,  j * 
J m=2  ’ J 


n n 2 n n 2 1/2 

xm-l,j)  + (ym,j  ‘ ym-l,j)  1 


(4.1.8) 


and 


ALi 


n+1 


n+1 

L2 


n 

L2 


n+1 

Ll 

ALj 

n 

Ll 


(4.1.9) 
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The  coordinates  yjj  appearing  in  those  FDEs  are  given  by  the 
following  equations: 


n+1  , n+1 

yij  - Lij  j-i 


(4.1.10) 


n+1 

yij 


AL-s 


,n+l 


n 


n+1 


ASij  + Lij 


ALi 


,n 


(4.1.11) 


j-2,3. . . ,JL 


where 


,n  JL 
ALi  - 2 

m=2 


n n 2 n n 2 1/2 

(xi,m  ' xi,m-l)  + (yi,m  " yi,m-l)  ] 


(4.1.12) 


,n  i n n 2 n 

ASij  "m^2  I <xi,m  - xi,m-l)  + (yi,m 


n 2 1/2 

yi.m-i)  ] 


(4.1.13) 


and 


,n+l  ,n+l 

, n+1  l2  * l1  , n 

ALj  ~ ' ALj 

,n  ,n 

l2  * 1*1 


(4.1.14) 
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n 


th 


In  the  above  equations , ALj  is  the  length  of  the  j grid  line 

th  n 

m the  physical  domain  at  the  n time  level;  ASjj  is  the  length 

_ , .th  . n+1 

of  the  j grid  line  from  i=l  to  i;  ALj  is  an  estimate  of  the 

th  th 

length  of  the  j grid  line  at  the  (n+1)  time  level.  The 

„ -n  ,n  , n+1 

physical  significance  of  AL* , ASjj  and  AL^  is  similar  to 

, n n • n+1 

that  explained  for  ALj , AS^j  and  ALj  . 

(9)  Apply  the  FDEs  derived  in  step  2 to  obtain  a "tentative"  solu- 
tion of  the  problem  being  solved  at  the  (n+1)  level  by  using 

, n n n+1  n+1 

the  grid  systems  given  by  (xjj  , yjj ) and  (Xij , ytj  ). 

(10)  Evaluate  the  smoothness  factor  by  using  the  grid  system  given  by 

n+1  n+1 

( xij . Yij  ) and  the  "tentative"  solution  obtained  in  step  9. 

This  is  very  involved  and  will  be  explained  in  Section  4.2. 

th 

(11)  Recalculate  the  locations  of  the  grid  points  at  (n+1)  time 

level  by  using  the  smoothness  factor  obtained  in  the  previous 

th 

step.  The  x-coordinates  of  the  grid  points  along  the  j grid 
th 


line  at  the  (n+1)  time  level  are  computed  by  using  Eqs. (3.3.5) 
to  (3.3.7).  In  these  equations,  xm=xmj ; the  smoothness  factor 
is  that  evaluated  at  step  10;  and  U(xm)-  U(xmj,ymj)  where  U is 

given  by  the  "tentative"  solution  obtained  in  step  9.  The  y- 

, . _ th 

coordinates  of  the  grid  points  along  the  i grid  line  are  ob- 
tained by  using  the  procedure  just  described  except  replace  x by 
y,  i by  j,  LX  by  L1(  L2  by  L2 , and  U(xmj,ymj)  by  U(xim,yim). 

(12)  Control  the  orthogonality  of  the  grid  system  obtained  in  step 
11.  This  is  very  involved  and  will  be  explained  in  Section  4.3. 

(13)  Recalculate  the  metric  coefficients  by  using  the  equations 

j . , . . n+1  n+1 

described  in  step  8 with  Xjj  and  y^j  obtained  in  step  12. 
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(14)  Use  the  FDEs  with  the  metric  coefficients  calculated  in  step  13 

. th 

to  obtain  the  "corrected"  solution  at  the  (n+1)  time  level. 

(15)  Replace  the  time  level  n by  n+1  and  calculate  the  total  time 
elapsed  by  using  Eq.  (3.3.11).  Repeat  steps  8 to  14  until  the 
total  time  elapsed  equals  to  the  duration  of  interest  specified 
in  step  5. 

This  completes  the  description  of  the  solution-adaptive,  grid 
generation  method  for  two-dimensional,  rectangularly  shaped  spatial 
domains  that  can  undergo  rectilinear  deformations.  In  the  next  two 
sections,  the  techniques  used  to  control  nonoverlapping  of  grid  lines 
and  the  smoothness  of  the  grid  system  as  well  as  the  procedure  used  to 
control  the  orthogonality  of  grid  system  are  explained. 

4,2  Nonoverlapping  of  Grid  Lines  and  Smoothness 

When  using  the  solution-adaptive,  grid  generation  method  presented 
in  the  last  section  to  generate  grid  points  inside  two-dimensional 
spatial  domains,  it  is  important  to  control  not  only  the  clustering  of 
grid  points  in  regions  where  most  needed,  but  also  the  nonoverlapping  of 
the  grid  lines  and  the  smoothness  of  the  grid  system  being  generated. 
In  this  investigation,  nonoverlapping  of  grid  lines  and  smoothness  were 
controlled  by  the  smoothness  factor  described  in  Section  3.2.2.  and  by 

monitoring  a cell  area  associated  with  each  grid  point. 

The  cell  area,  Ag,  associated  with  each  grid  point  can  be  defined  by 

the  magnitude  of  the  cross  product  of  the  vectors  tangent  to  the  curves 

£ and  rj  at  grid  point  P shown  in  Fig  4.2  and  represented  by  the  follow- 
ing equation: 
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Figure  4.2  Definition  of  the  vectors  tangent  to  the  curves  £ 
and  r)  at  grid  point  P. 
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Ag 


dr  dr 

A r?  x Af 

dr,  dt 


Ag 


-*•  3y  3x  _►  dy  .+ 

Ar?i  + Ar?j)  x ( Afi  + A£j ) 

drj  dr/  d{  d£ 


(A. 2.1) 


which  after  efetuating  the  vectorial  operations  reduces  to 

dx  dy  dx  dy 

Ag  (4.2.2) 

dr)  dr)  d£ 

-4 

when  Ar)  — A£  — 1 . In  Eq.  (4.2.1)  i and  j are  unit  vectors  point  in  the 
x-  and  y-direction,  repectively. 

The  equation  above  is  the  Jacobian  of  the  coordinate  transformation 
given  by  Eq.  (2.2)  when  t=r . This  fact  will  be  useful  because  the 
Jacobian  has  already  been  calculated  while  computing  the  "tentative" 
solution.  The  Jacobian  can  be  evaluated  numerically  by  either  one  of 
the  following  two  ways : 

Agij  - Jij  ~ t(xi+lj  - xi-l , j ) (yi ; j+i  - yij-l)  - 

(xi,j+l  - xi,j-l)(yi+l,j  - yi-l,j)]/4  (4.2.3) 

or 

Agij  “ Jij  * t (xi+l/2 , j - xi-l/2 , j ) (yi , j+1/2  - yi.j-1/2)  * 

<xi,j+l/2  ' xi , j - l/2> (yi+1/2 , j - yi-l/2,j>]  (4.2.4) 

The  cell  areas  represented  by  these  two  equations  are  depicted  in 
Fig.  4.3.  In  this  investigation,  Eq.  (4.2.3)  is  used. 
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(b) 


Figure  4.3  Correspondence  between  the  Jacobian  and  the  cell  area. 

(a)  The  shaded  area  is  the  cell  area  defined  by  Eq.  (4.2.3) 

(b)  The  shaded  area  is  the  cell  area  defined  by  Eq.  (4.2.4) 
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With  the  concept  of  cell  area  defined,  the  algorithm  for  controlling 
smoothness  can  be  described  in  the  following  steps: 

(1)  Set  j=2. 

n n th 

(2)  Determine  Jmax  and  Jm£n  along  the  j grid  line. 

n n th 

(3)  Calculate  wmax  and  wm£n  along  the  j grid  line  by  using  Eq. 

th 

(3.2.14)  where  U is  the  solution  at  the  n time  level. 
n+1  n+1  th 

(4)  Calculate  wmax  and  wmin  along  the  j grid  line  by  using  Eq. 

ttl 

(3.2.14)  where  U is  the  tentative  solution  at  (n+1)  time  level. 

n+1 

(5)  Calculate  the  ratio  (Jmax/Jmin)  bY  using  Eq.  (3.2.16)  where 
Amax/Amin  are  replaced  by  Jmax/Jmin- 

(6)  Calculate  the  smoothness  factor  SF  by  using  Eq.  (3.2.13)  with 

n+1  n+1 

Amax  anb  Amin  replaced  by  Jmax  and  Jmin,  respectively. 

(7)  Replace  j by  j+1  and  repeat  steps  2 to  6 until  j-JL-1. 

th 

(8)  Repeat  step  1 to  7 for  the  i grid  line  from  i-2  to  i-  IL-1. 

n+1 

(9)  Calculate  the  Jjj  . 

The  grid  system  generated  by  the  following  steps  1-9  described  above 
could  contain  overlapping  grid  lines.  The  sufficient  condition  for  non- 
overlapping grid  lines  in  the  physical  domain  is  that  the  cell  area 

associated  with  each  grid  point  (or  the  Jacobian)  never  vanishes. 

T,  . , . . , n+1  n+1 

ihis  can  be  achieved  by  controlling  the  value  of  Jjj  . If  J ± j calcu- 
lated in  step  9 is  greater  than  Jcrt  (a  value  given  a priori) , then  the 

th 

algorithm  can  be  stopped.  If  J^j  calculated  in  step  9 is  less  than 
^crt>  then  the  algorithm  continues  from  step  9 with  the  following  step: 

(10)  Rewrite  Eq.  (3.2.16)  as  follows: 
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n+1 


“max 


^max 

^max 

wmin 

^min 

^min 

wmax 

n 


n 


n+1 


(4.2.5) 


wmin 


(11)  Repeat  steps  2-6  for  all  i and  j grid  lines  where  Jij<Jcrt 
with  one  exception.  The  exception  is  in  step  5 instead  of  using 
Eq.  (3.2.16),  use  Eq.  (4.2.5). 

This  completes  the  algorithm  for  controlling  nonoverlapping  of  grid 
lines  and  smoothness.  In  the  next  section  a technique  for  controlling 
the  orthogonality  of  the  grid  system  is  presented. 


4.3  Orthogonality  of  the  Grid  System 

As  noted  earlier,  grid  systems  for  two-dimensional  spatial  domains 
should  be  nearly  orthogonal.  In  Chapter  I,  it  was  explained  that  nearly 
orthogonal  means  that  the  grid  lines  should  intercept  each  other  with  an 
angle  approximately  between  45  and  135  degrees.  This  fact  can  be  used 
to  control  near  orthogonality  as  shown  in  Fig  4.4. 

In  Fig.  4.4,  the  locations  of  grid  points  A and  B have  already  been 
determined  and  the  lines  Aa,  Ab,  Be  and  Bd  define  a sector  of  the 
physical  domain  in  which  the  coordinates  of  the  new  grid  point  to  be 
generated  can  be  located  to  satisfy  the  minimum  requirements  for  the 
grid  to  be  nearly  orthogonal. 

The  lines  la  and  lb  emanating  from  point  A can  be  represented  by  the 


following  equations: 
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Figure  4.4  Orthogonality  of  grid  lines  in  the  physical  domain. 


54 


Line  la 


y - 

mla  x + hla 

(4.3.1) 

mla  “ 

tan(  arctan  m^  + ir/4  ) 

(4.3.2) 

hla  “ 

yi  - mla  X1 

(4.3.3) 

Line  lb 


y 

“ mlb  x + hlb 

(4.3.4) 

mlb 

= tan(  arctan  m^  + 3tt/4  ) 

(4.3.5) 

hlb 

- yi  - mib  xi 

(4.3.6) 

where  mb  is  the  tangent  to  the  grid  line  at  grid  point  A. 

Similarly,  lines  2c  and  2d  emanating  from  grid  point  B can  be  re- 
presented by  the  following  equations: 

Line  2c 

y “ m2c  x + h2c  (4.3.7) 

m2c  “ tan(  arctan  m2  + tt/4  ) (4.3.8) 

h2c  - y2  * m2c  x2  (4.3.9) 


Line  2d 


y - m2d  x + h2d 


(4.3.10) 


m2d  “ tan(  arctan  m2  + 3?r/4  ) (4.3.11) 

h2d  “ Y2  ' 


m2d  x2 


(4.3.12) 
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where  m2  is  the  tangent  to  the  grid  line  at  grid  point  B. 

Finally,  to  accomplish  the  requirements  of  near  orthogonality,  the 
coordinates  of  the  grid  point  to  be  generated  have  to  satisfy  the 
following  inequalities: 

y > mla  x + hla  (4.3.13) 
y > m2c  x + h2c  (4.3.14) 
x > (y  * hlb>/mlb  (4.3.15) 
x > (y  - h2d)/m2d  (4.3.16) 

The  inequalities  expressed  by  Eqs . (4.3.13)  to  (4.3.16)  were  used 
to  locate  the  coordinates  of  the  grid  points  in  the  physical  domain  to 
ensure  the  grid  generated  would  be  nearly  orthogonal  during  the  appli- 
cation of  the  two-dimensional  grid  generation  method  developed. 

This  concludes  the  description  of  the  solution- adaptive  grid  genera- 
tion method  for  two-dimensional,  rectangularly  shaped  spatial  domains. 
Several  examples  illustrating  how  the  method  developed  in  this  chapter 
can  be  applied  for  generating  two-dimensional  grid  systems  are  presented 
in  Chapter  VIII . 


CHAPTER  V 


THE  SAAGG  METHOD  FOR  2-D,  ARBITRARILY  SHAPED  DOMAINS 

The  two-dimensional,  solution- adaptive , algebraic  grid  generation 
(SAAGG)  method  developed  in  Chapter  IV  can  be  used  with  any  non- 
adaptive,  grid  generation  technique  to  generate  solution- adaptive  grid 
systems  in  arbitrarily  shaped  2-D  spatial  domains  that  can  deform  in  an 
arbitrary  manner.  This  is  demonstrated  in  this  investigation  by  using 
the  method  developed  in  Chapter  IV  in  conjunction  with  the  Two -Boundary 
Technique.  This  demonstration  is  presented  in  the  following  order. 
First,  the  Two-Boundary  Technique  is  briefly  described.  Afterwards,  the 
criteria  for  controlling  the  nonoverlapping  of  grid  lines  in  the  physi- 
cal domain  and  the  orthogonality  of  the  grid  system  at  the  boundaries 
are  discussed.  Finally,  how  the  method  developed  in  Chapter  IV  is  in- 
corporated into  the  Two -Boundary  Technique  is  described. 

5,1  Overview  of  the  Two -Boundary  Technique 
The  Two -Boundary  Technique  is  a non- solution- adaptive , algebraic 
grid  generation  method  based  on  interpolation  techniques  [26].  In  gen- 
eral, Hermite  interpolation  is  applied  between  two  boundaries  of  the 
spatial  domain  to  establish  the  coordinate  transformation  which  maps  the 
physical  domain  onto  the  computational  domain  [27].  Here,  it  is  impor- 
tant to  note  that  the  Two -Boundary  Technique  only  maps  two  boundaries  of 
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the  spatial  domain  correctly.  The  other  two  boundaries  of  the  spatial 
domain  can  be  mapped  correctly  only  if  they  are  straight  lines  [3,18]. 
In  practice,  there  are  some  problems  in  which  it  is  necessary  to  map 
four  boundaries  of  the  spatial  domain  correctly.  In  these  cases,  the 
Four-Boundary  Technique  can  be  used  instead  of  the  Two-Boundary 
Technique  [3], 

The  Two -Boundary  Technique  can  be  explained  by  considering  Fig. 
(5.1)  where  the  correspondence  between  the  boundaries  of  the  physical 
and  the  computational  domains  is  shown.  In  that  figure,  the  two 
boundaries  AD  and  BC  of  the  spatial  domain  are  described  by  the 
following  parametric  equations: 

X1  “ xl(f.»?-0,r)  yi  - y]_(£,»?-0,r)  (5. 1.1, 2) 

x2  “ x2(£>’7“1>t)  yi  ~ y2 (£,»?-l,r)  (5. 1.3, 4) 

where  the  subscripts  1 and  2 denote  boundaries  AD  and  BC,  respectively. 
Note  that  since  x^ , X2 , y]_ , y2  are  functions  of  r,  the  boundaries  AD 
and  BC  can  deform  with  time. 

The  Two-Boundary  Technique  based  on  Hermite  interpolation  gives  the 
following  coordinate  transformation  between  the  physical  domain  in  x-y-t 
coordinate  system  and  the  computational  domains  in  coordinate 
system  [3,28] : 


x(£. *J.O  - X]_(£,r)hi(»j)  + ,r)h2(v)  + 


5x(£,t7=0  ,r) 


dx(£,»7-l,r) 


dr) 


dr) 


■*13  (r?)  + 


'^4  (»7) 


(5.1.5) 
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y 


(a) 


(b) 


Figure  5.1  Correspondence  between  the  boundaries  of  the  physical 
domain  and  the  computational  domain. 

(a)  Arbitrarily  shaped  physical  domain  (x-y-t). 

(b)  Computational  domain  (£-r?-r)  . 
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y(£,»?,0  - yi(Z  ,r)hi(r))  + y2(£,r)h2  (r?)  + 

3y(^,»?-0,T)  ay(^,»7-l,T) 

h3(T7)  + h4(rj)  (5.1.6) 

dr]  dr] 

The  coefficients  h^,  h2 , I13 , and  h4  in  the  equations  above  are  given 
by  the  following  expressions  [3]: 


hi(»?)  - 2»73  - 3>72  + 1 (5.1.7) 

h2(r?) 2»73  + 3 r?2  (5.1.8) 

h3 (»7)  - r?3  - 2»j2  + r?  (5.1.9) 

h4(r?)  “ *?3  - *?2  (5.1.10) 


The  partial  derivative  terms  appearing  in  Eqs . (5.1.5)  and  (5.1.6) 
are  chosen  so  that  grid  lines  intersect  the  boundaries  AD  and  BC 
perpendicularly  in  the  physical  domain.  Accordingly,  these  partial 
derivative  terms  are  given  as  follows  [3]: 

3x(£,r7«0,r)  dyi 

- - KX(0  (5.1.11) 

dt]  5? 


9y(C.»7“0 ,0 

dr] 


3xi 


“ + KX(0 


(5.1.12) 


3x(£,r/=l,r) 


(5.1.13) 


dy2 

- - k2(0  

9ri  3| 


3y(C.»7“l  ,t) 

dr) 


3x2 

+ k2(0  * 

a? 


(5.1.14) 


The  coefficients  and  K2  appearing  in  the  equations  above  are 
determined  by  trial  and  error  to  ensure  nonoverlapping  of  the  grid  lines 
in  the  physical  domain.  In  order  to  incorporate  the  solution-adaptive 
method  developed  in  Chapter  IV  into  the  Two-Boundary  Technique,  these 
coefficients  must  be  determined  analytically. 

5,2  Extended  Two-Boundary  Technique 
In  this  section,  the  Two -Boundary  Technique  is  extended  so  that  it 
can  be  applied  with  the  method  developed  in  Chapter  IV.  First,  a 
procedure  is  described  for  expressing  the  coefficients  K]_  and  K2  in 
analytical  form.  Second,  criteria  are  defined  for  ensuring  that  grid 
lines  do  not  overlap  each  other  or  the  boundaries.  Third,  a technique 
is  presented  which  ensures  that  the  angles  at  which  grid  lines  intersect 
the  two  boundaries  are  within  certain  limits. 

5,2.1.  Analytical  Expressions  for  K1  and  K? 

The  presentation  of  the  procedure  for  determining  the  coefficients 
K1  and  k2  is  facilitated  by  substituting  Eqs . (5.1.11)  to  (5.1.14)  into 
Eqs . (5.1.5)  and  (5.1.6)  as  shown  below. 
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x(£,»?,r)  - 

xlhi 

+ x2h2 

- KxBxh3  - K2B2h4 

(5.2.1) 

y (£>*?> r ) - 

yihi 

+ y2h2 

+ KxAxh3  + K2A2h4 

(5.2.2) 

where  the  partial  derivatives 

appearing  in  Eqs.  (5.1.11) 

to  (5.1.14) 

were  represented  by 

3xx 

dx2 

AX  - - 

a2  - 

(5. 2. 3, 4) 

a£ 

ac 

3yi 

ay2 

Bx  - - 

b2  - 

(5. 2. 5, 6) 

5? 

a? 

In  order  to  determine  analytical  expressions  for  Kx  and  K2 , first,  a 
grid  line  should  be  specified,  a priori,  inside  the  physical  domain.  In 
this  investigation,  the  locations  of  grid  points  along  this  grid  line 
were  chosen  to  be 


xi  + x2 

x(£,r?-0.5)  - Q(0  - € + (5.2.7) 

2 


yi  + yi 

y(£,»J-0.5)  « PCO  - yi  + (5.2.8) 

2 


where  X]_,  x2 , yx,  and  y2  are  functions  of  £ and  are  given  by  Eqs . 
(5.1.1)  to  (5.1.4),  and  e is  a small  quantity  between  ± O.lAx. 
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Next,  the  values  of  h^,h2,  h3,  and  h4  needed  to  be  determined. 
Since  r?- 0.5,  Eqs . (5.1.7)  to  (5.1.10)  give 

1 1 

hi  - h2  - and  h3  - -h4  = (5.2.9,10) 

8 8 

Substituting  Eqs.  (5.2.7)  to  (5.2.10)  into  Eqs.  (5.2.1)  and  (5.2.2), 
and  solving  for  the  coefficients  and  K2  yield  the  following  expres- 
sions : 


Kl 


8 (QA2  + PB2)  - 4[A2(x^  + x2)  + B2 (yj_  + y2)] 
(A3B2  * A2Bl) 


(5.2.11) 


8 (QAq  + PB^)  - 4[A]_(xx  + x2)  + Bl(yi  + Y2)] 

K2  - (5.2.12) 

(A1B2  - A2Bq ) 

The  equations  above  are  the  desired  analytical  expressions  for  the  coef- 
ficients and  K2. 

Here,  it  is  noted  that  the  expressions  for  the  coefficients  and 
K2  given  by  Eqs.  (5.2.11)  and  (5.2.12)  approach  infinity  when  the  deno- 
minator in  those  equations  approaches  zero.  This  problem  can  be  over- 
come by  setting  and  K2  equal  to  the  average  values  of  the  coef- 

ficients immediately  before  and  after  the  grid  point  in  which  Eqs. 
(5.2.11)  and  (5.2.12)  become  singular. 
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The  analytical  expressions  for  Kx  and  K2  defined  by  Eqs . (5.2.11) 
and  (5.2.12)  will  be  useful  when  the  solution- adaptive  grid  generation 
method  developed  in  Chapter  IV  is  used  with  the  Two-Boundary  Technique. 
The  influence  of  the  coefficients  Kx  and  K2  in  controlling  the  distri- 
bution of  the  grid  points  in  the  physical  domain  is  studied  in  the  next 
section.  This  study  is  made  to  establish  the  criteria  for  nonoverlap- 
ping grid  lines. 

5.2.2  Criteria  for  Nonoverlapping  Grid  Lines 

The  grid  r;-grid  lines  in  the  physical  domain  are  determined  by  using 
Eqs.  (5.2.1)  and  (5.2.2)  and  setting  r?  to  a constant  between  0 and  1 in 
Eqs.  (5.1.7)  to  (5.1.10).  The  values  of  Kx  and  K2  determine  whether  any 
given  r?-grid  line  will  be  located  closer  or  farther  away  from  the 
boundaries.  Thus,  it  is  important  to  determine  the  limits  of  these 
coefficients  in  order  to  control  the  nonoverlapping  of  grid  lines. 

The  r/-grid  lines  will  not  overlap  the  two  boundaries  AD  and  BC  of 
the  physical  domain  if  the  y-coordinates  of  those  grid  lines  satisfy  the 
following  inequalities: 


By  using  Eqs.  (5.1.6)  to  (5.1.9)  and  setting  rj— 1— Arj , the  following 
equations  can  be  written: 


y(£,Ar?)  > yi(£,0) 


(5.2.13) 


y(*,l-A.j)  < y2(£ , 1) 


(5.2.14) 


3 


2 


hX  - 2(1-Ar?)  - 3(1-Aij)  + 1 


(5.2.15) 
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3 2 

h2  - -2(1-Ary)  + 3(1-Ar;) 

(5.2.16) 

3 2 

I13  - (1-Ary)  - 2(1-Ary)  + (1-Ary) 

(5.2.17) 

3 2 

h4  -=  (1-Ary)  - (1-Ary) 

(5.2.18) 

By  neglecting  the  higher  order  terms  in  Ary , the  equations  above  give 
hq=0,  h2=l,  h3«0,  and  I14 — A rj . By  substituting  these  results  into  Eq. 
(5.2.2)  and  the  resultant  equation  into  Eq.  (5.2.14),  the  following 


inequality  can  be  written: 

K2  Ary  A2  > 0 

(5.2.15) 

Since  Ary  is  always  positive  and  the  partial  derivative  A2  is  also 
positive,  the  only  solution  which  satisfies  the  inequality  above  is 


K2  > 0 

Similarly,  for  ry=Ary  the  following  inequality  can  be  found: 

(5.2.16) 

0 

A 

rH 

(5.2.17) 

Thus,  the  criterion  for  the  grid  lines  not  to  overlap  the  boundaries 
AD  and  BC  of  the  physical  domain  is  that  the  coefficients  K]^  and  K2  are 
always  positive. 
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The  criteria  for  nonoverlapping  r?-grid  lines  in  the  interior  of  the 
physical  domain  are  very  involved  and  are  described  below  for  spatial 
domains  which  contains  a line  of  symmetry.  For  those  spatial  domains, 
it  can  be  deduced  from  Eqs . (5.2.1)  to  (5.2.4)  that  A1-A2-A,  and 
Bi— B2-B.  Similarly  from  Eqs.  (5.2.11)  and  (5.2.12),  it  can  be  deduced 
that  Ki=K2“K. 

To  establish  the  criteria  for  nonoverlapping  of  interior  r/-grid 
lines,  it  is  necessary  to  investigate  the  influence  of  the  coefficient  K 
on  the  distribution  of  grid  points.  This  influence  was  analyzed,  in 
this  investigation,  through  the  ratio  between  Am^  and  Ab^  determined  by 
the  coordinates  of  grid  points  at  r,  s,  t and  m,  along  the  grid  line  £ ^ , 
shown  in  Fig.  5.2.  This  ratio  is  given  by 


Am^ 

Ab^ 


(5.2.18) 


where 


Amj;  -= 


[ (xm 


? 0 V2 

xt)  + (ym  - yt)  k 


(5.2.19) 


and 


Abi  - ± [(xs  - xr)2  + (ys  - yr) 2 ] J/2 


(5.2.20) 


In  the  above  equations,  the  negative  sign  is  used  to  indicate  those 
cases  in  which  (ym  - yt)  <0  or  (ys  - yr)  < 0.  The  coordinates  of 
the  grid  points  at  s,  t and  m along  grid  line  are  given  by  Eqs. 
(5.2.1)  and  (5.2.2)  by  setting  r?“=Ar/  , r?=0.5-Ar7,  and  tj— 0.5,  respectively. 
The  coordinates  of  the  grid  point  at  r along  the  boundary  AD  are  given 
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V = .5 


0 x 


Figure  5.2  Grid  lines  in  a spatial  domain  which 
contains  a line  of  symmetry. 
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by  Eqs . (5.1.1)  and  (5.1.2).  The  substitution  of  the  coordinates  thus 
obtained  into  Eqs.  (5.2.19)  and  (5.2.20)  gives  an  expression  for  Am^  and 
Ab^  as  function  of  and  . Hence,  the  ratio  oc^  can  be  expressed  in 
terms  of  and  as  follows: 

ai  “ a(£i,Ki)  (5.2.21), 

The  exact  functional  form  of  Eq.  (5.2.21)  is  not  shown  because  it 
depends  on  the  parametric  equation  used  to  describe  the  boundary  AD  of 
the  spatial  domain.  Figure  5.3  shows  a typical  relationship  between  oc 
and  K for  several  grid  lines  In  Fig.  5.3,  also  three  regions  can  be 
distinguished.  When  oc  > 1 , ry- grid  lines  tend  to  cluster  near  the 
boundaries.  When  0 < a < 1,  ry- grid  lines  tend  to  stretch  away  from 
the  boundaries.  Finally,  when  cc  < 0 , ry-  grid  lines  overlap. 

A few  other  observations  can  also  be  made  from  Fig.  5.3.  First,  if 
K < Ka,  then  every  grid  point  along  ry-grid  lines  in  the  physical  domain 
is  clustered  to  the  boundaries.  Second,  if  Ka  < K < Kb , then  some  grid 
points  along  ry-grid  lines  will  be  clustered  near  the  boundary  while 
others  will  be  stretched  away  from  the  boundary  depending  upon  the  value 
of  oc^.  Third,  if  Kb  < K < Kc,  then  every  grid  point  along  ry-grid  line 
will  be  stretched  away  from  the  boundary.  Finally,  if  K > Kc , then  some 
ry-grid  lines  will  be  overlaped. 

The  limits  Ka,  Kb,  and  Kc  should  be  determined  before  the  solution- 
adaptive  grid  generation  method  developed  in  Chapter  IV  is  to  be  applied 
with  the  Two-Boundary  Technique.  This  is  because  the  locations  of  the 
grid  points  in  physical  domain  are  influenced  by  the  coefficients  K as 
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0 1 ai 


Figure  5.3  Relationship  between  the  ratio  ex^-Am^/Abj. 
and  the  coefficient  K. 
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well  as  by  the  stretching  function  defined  in  Chapter  II. 

Thus,  the  criterion  for  nonoverlapping  grid  lines  is  that  the  ratio 
ai  be  positive  for  all  of  the  grid  lines  in  the  physical  domain;  i.e., 

Ami 

«i  - > 0 (5.2.22) 

At>i 

It  was  shown  in  this  section  that  the  coefficient  K must  be  positive 
for  the  grid  lines  to  not  overlap  the  two  boundaries  AD  and  BC  of  the 
physical  domain.  Thus,  if  K is  positive,  then  Ab^  must  be  positive 
which  implies  from  Eq.  (5.2.22),  that  Am^  must  be  also  positive. 

From  Eq . (5.2.19),  it  can  be  seen  that  the  condition  for  Am^  to  be 

positive  is  that 

ym  - Yt  > 0 (5.2.23) 

Since  f?- 0.5  is  a line  of  symmetry,  ym  in  the  equation  above  is 
determined  from  the  following  expression; 

D 

ym  - yi  + — (5.2.24) 

where  y^  is  obtained  from  Eq.  (5.1.2)  and  D/2  is  the  distance  shown  in 
Fig.  5.2.  In  Eq.  (5.2.23),  yt  is  obtained  from  Eq.  (5.2.2).  In  this 
equation  the  coefficients  h^,  h2 , 113  and  h^  are  determined  from  Eqs . 
(5.1.7)  to  (5.1.10)  by  setting  r;=0.5-Ar;  and  neglecting  higher  order 
terms  in  A rj . After  simplification,  the  following  equations  are  found: 
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3 1 

h]_  - A t]  + (c 

2 2 

3 1 

h2  " - Ary  + (5 

2 2 

1 A ry 

h3  " + (5 

8 4 

1 A r) 

h4 +.  (5 

8 4 


The  substitution  of  Eqs . (5.2.25)  to  (5.2.28)  into  Eq.  (5.2.2) 


, 3 1 1 
yt  “ yi(— A r,  + — ) + (yi  + D)(— 


3 A r) 

—Ary)  + (KA) 


D A»y 

yt  - yi  + - (3D  . ka) 


(5 


Finally,  the  substitution  of  Eqs.  (5.2.24)  and  (5.2.29)  into  Ec 
(5.2.23)  gives  the  following  expression: 


Ary 

(3D-KA)  > 0 

2 


(5. 


Since  A »y  is  positive  and  K is  also  positive,  the  criterion  for 
overlapping  ry-grid  lines  in  the  physical  domain  becomes 


.2.25) 


.2.26) 


.2.27) 


.2.28) 


gives 


.2.29) 


2.30) 


non- 
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3D 


0 < K < 


(5.2.31) 


A 


Here,  it  is  noted  that  the  inequality  K > 0 is  valid  for  all  types 
of  spatial  domains  but  the  second  inequality  K < 3D/A  is  only  valid  for 
those  spatial  domains  which  possess  one  line  of  symmetry  between  two 
boundaries  AD  and  BC. 

5-2.3 Criteria  for  Orthogonality  at  Boundaries 

By  using  the  Two-Boundary  Technique  with  Hermite  interpolation,  Eqs . 

(5.1.11)  to  (5.1.14)  ensure  that  the  grid  lines  intersect  orthogonally 

the  boundaries  of  the  physical  domain.  However,  this  is  not  true  after 

the  continuous  grid  lines  have  been  replaced  by  a finite  number  of  grid 

points.  In  this  section,  the  Two-Boundary  Technique  is  extended  to 

ensure  that  the  deviation  from  orthogonality  at  the  boundaries  of  the 

physical  domain  is  within  certain  limits. 

The  deviation  from  orthogonality  at  the  boundaries  (0)  is  shown  in 

Fig.  5.4.  In  this  figure  it  can  be  seen  that  the  angle  0 is  formed  by 

-*■  -*• 

the  vector  normal  to  the  boundary  e^  and  the  vector  a. 

The  vector  normal  to  the  boundary  can  be  expressed  by 


3x(Ci,r/-0) 


3y(?i,r?-0)  _ 


i + 

a? 


j 

dt 


(5.2.32) 


where  the  x and  y are  given  by  Eqs.  (5.2.1)  and  (5.2.2),  respectively. 
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Figure  5.4  Deviation  from  orthogonality  at  the  boundaries 
of  the  physical  domain. 


The  vector  a is  expressed  in  terms  of  the  coordinates  of  grid  points 


P(x,y)  and  P0(x0,y0)  shown  in  Fig.  5.4,  as  follows: 


a 


(x  - xQ)  i + (y  - yQ)  j 


(5.2.33) 


The  coordinates  of  grid  point  P(x,y)  are  determined  by  Eqs.  (5.2.1) 
and  (5.2.2)  for  a given  grid  line  £i  and  with  rj-Ar?  as  function  of  the 
coefficients  K^Ci)  and  K2(^). 

The  coordinates  of  grid  point  P0(xo>yo)  at  the  boundary  of  the 
physical  domain  are  determined  by  Eqs . (5.1.1)  and  (5.1.2)  once  the  grid 
line  has  been  specified. 

The  angle  8 between  the  vectors  e^  and  a is  obtained  from  analytical 
geometry  and  is  given  by 


By  substituting  Eqs.  (5.2.32)  and  (5.2.33)  into  the  equation  above  and 
effectuating  the  vectorial  operation,  the  following  equation  is 


8 - cos 


-1  e£  . a 


(5.2.34) 


obtained: 


-1 


ElCi  + 


cos 


(5.2.35) 


2 2 2 2 1/2 
[(Ei  +FX  )(Ci  +Di  )] 


where  the  following  substitutions  were  made: 
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El 


dx(ti>v-0) 


3Z 


El 


dy(Zi.r)-O) 


dt 


(5.2.36,37) 


C1  “ <x  * xo)  &1  - (y  - yQ) 


(5.2.38,39) 


Equation  (5.2.35)  is  a function  of  £it  Ar? , and  the  coefficients 
and  K2 ; i . e . , 


*i  - tf«i,A^,K1,K2) 


(5.2.40) 


The  exact  functional  form  of  Eq.  (5.2.40)  is  not  shown  because  it 
depends  on  the  parametric  equation  used  to  describe  the  boundaries  of 
the  spatial  domain.  Figure  5.5  shows  a typical  relationship  between  6 
and  for  a given  value  of  r?=Ar/  and  several  values  of  K. 

Figure  5.5  shows  that  the  bigger  the  coefficient  K the  smaller  is 
the  deviation  from  orthogonality  at  the  boundaries.  Unfortunately, 
constraints  imposed  by  Eq.  (5.2.31),  to  ensure  the  grid  lines  do  not 
overlap,  sometimes,  do  not  allow  the  coefficient  K to  be  large  enough 
to  keep  the  deviation  from  orthogonality  within  the  limits  required  for 
implementation  of  the  boundary  conditions. 

In  this  investigation,  two  techniques  were  developed  for  controlling 
the  orthogonality  at  the  boundaries.  The  first  technique  controls  or- 
thogonality by  specifying  the  minimum  number  of  grid  points  JL,  in  y- 
direction,  needed  to  keep  the  deviation  from  orthogonality  within  a 
prescribed  limit.  The  second  technique  controls  the  orthogonality  by 
determining  a stretching  function  which  redistributes  existing  grid 
points  to  keep  the  deviation  within  the 


prescribed  limit. 
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Figure  5.5  Deviation  from  orthogonality  at  the  boundaries 
of  the  physical  domain  as  function  of  and 

the  coefficient  K. 
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Figure  5.6  shows  how  the  deviation  from  orthogonality  can  be  reduced 
by  clustering  the  grid  lines  to  the  boundaries  of  the  physical  domain. 
In  this  figure,  the  angle  6 is  the  deviation  from  orthogonality  at  the 
boundary  before  stretching  and  p is  the  maximum  deviation  from  ortho- 
gonality allowed. 

To  reduce  the  deviation  from  orthogonality  by  the  first  technique, 
Eq.  (5.2.35)  is  rewritten  to  introduce  the  maximum  deviation  from 
orthogonality  allowed,  as  shown  below: 

-1  Elcl  + F]_Di 

0 “ cos  (5.2.41) 

2 2 2 2 1/2 
[(E1  +F1  )(C1  +DX  )] 

Once  the  maximum  coefficient  K has  been  determined  by  Eq.  (5.2.31) 
everything  on  the  right  hand  side  of  Eq.  (5.2.40)  is  known  except  for 
the  variables  and  D]_.  These  variables  are  defined  by  Eqs . (5.2.36) 
and  (5.2.37)  which  are  expressed  as  functions  of  the  and  A r) . Thus, 
for  the  maximum  angle  of  deviation  p and  a given  grid  line  a new 
value  of  Ar?,  denoted  by  Ar?*,  can  determined  from  Eq.  (5.2.41). 

The  grid  spacings  Ar?  and  Ar?*  are  related  to  the  number  of  grid  lines 
by  the  following  expressions: 

1 1 

Ar?  - Ar?*  . (5.2.42,43) 

JL  - 1 JL*  - 1 

Now,  the  minimum  number  of  grid  points  which  will  keep  the  deviation 
from  orthogonality  no  greater  than  p is  determined  by  substituting  the 
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Figure  5.6  Reduction  of  deviation  from  orthogonality 
by  using  stretching  function. 
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value  of  Art*  obtained  from  Eq.  (5.3.41)  in  Eq.  (5.3.43),  as  follows: 


1 

JL*  - 1 + (5.2.44) 

Ar,* 

The  second  technique  to  reduce  the  deviation  from  orthogonality  was 
developed  for  problems  in  which  the  number  of  grid  points  cannot  be 
changed.  For  these  problems,  a stretching  function  must  be  determined 
in  order  to  redistribute  the  grid  points  in  the  physical  domain  and  move 
the  grid  lines  toward  the  boundaries. 

To  demonstrate  the  second  technique,  let  the  stretching  function  be 
a second-order  polynomial  shown  below: 


* 2 

V - At?  + Br?  + C (5.2.45) 


where  r?  and  rj*  are  the  r?- coordinates  before  and  after  stretching.  The 
coefficients  A,  B,  and  C in  Eq.  (5.2.45)  are  defined  by: 


(JL  - 1)  - 2 (JL  - 1) 

B - 

(G  - 3)  (JL*-1) 


(5.2.46) 


A - 2(1  - B) 


(5.2.47) 


C = 0 


(5.2.48) 


when  r?  < 0.5  and 


2 * 

2 (G  - 1)  - 3 (G  - 1) 


B - 


(G  - 1)  (2G  - 1) 


(5.2.49) 
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2 

A - (1  - B) 

3 

1 

C - (1  - B) 

3 

when  r,  > 0.5.  In  Eqs.  (5.2.46)  through  (5.2.51) 
tutions  were  made: 


(5.2.50) 


(5.2.51) 
the  following  substi- 


JL  - 2 

G = 

JL  - 1 


and 


(5.2.52,53) 


The  stretching  function  defined  by  Eq.  (5.2.45)  redistribute  the 
grid  points  in  the  physical  domain  so  that  deviation  from  orthogonality 
was  no  greater  than  p. 

This  concludes  the  description  of  the  extended  Two -Boundary 
Technique.  In  the  next  section,  an  algorithm  for  applying  the  solution- 
adaptive  grid  generation  method  developed  in  Chapter  IV  with  the  ex- 
tended Two-Boundary  Technique  is  presented. 


— Algorithm  for  2-D,  Solution-Adanti ve  Grid  Systems 
This  section  presents  the  algorithm  for  generating  solution-adaptive 
grid  systems  for  two-dimensional  problems  in  which  the  boundaries  of  the 
physical  domain  can  move  or  be  stationary.  The  algorithm  will  map  the 
grid  points  between  the  physical  domain  in  x-y-t  coordinates  system  and 
the  computational  domain  in  £-»?-r  coordinates  system  by  using  the  SAAGG 
method  present  in  Chapter  IV  in  conjunction  with  the  Two -Boundary 
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Technique.  The  algorithm  is  as  follows: 

(1)  Transform  the  governing  equations  from  the  coordinate  system  of 
the  physical  domain  to  the  coordinate  system  of  the  computa- 
tional domain. 

(2)  Apply  any  suitable  finite-difference  method  to  express  the 
transformed  governing  equations  as  a set  of  finite-difference 
equations  (FDEs) . 

(3)  Specify  the  number  of  grid  lines  IL  in  the  ^-direction  and  the 
number  of  grid  line  JL  in  the  r\ -direction. 

(4)  Set  the  initial  time  level,  n,  to  zero  (which  corresponds  to 
time  t— 0)  and  specify  the  duration  of  interest  T (i.e.  0<t<T) . 

(5)  Decide  the  correspondence  between  the  boundaries  of  the  physical 
domain  and  the  boundaries  of  the  computational  domain. 

(6)  Express  the  coordinates  of  two  boundaries  of  physical  domain  in 
terms  of  the  boundary- fitted  coordinates  as  shown  by  Eqs. (5.1.1) 
to  (5.1.4). 

(7)  Determine  the  partial  derivatives  appearing  in  Eqs.  (5.2.1)  and 
(5.2.2)  by  using  expressions  obtained  in  step  6. 

(8)  Specify  the  locations  of  the  grid  points  in  the  computational 
domain. 


i - 1 

£ i “ 

IL  - 1 


j - 1 

V±  - 


i-1 , 2 , . . . , IL  (5.2.54) 


JL  - 1 


j-1,2 JL 


(5.2.55) 
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(9)  Calculate  the  coefficients  h^ , ^2’  ^3  and  by  using  Eqs . 

(5.1.7)  to  (5.1.10). 

(10)  Calculate  the  coordinates  of  grid  points  given  by  the  functions 
P and  Q shown  in  Eqs.  (5.2.7)  and  (5.2.8). 

(11)  Determine  the  coefficient  K by  using  P and  Q calculated  in  step 
10  and  Eqs. (5.2.11)  and  (5.2.12). 

(12)  Use  the  coefficients  calculated  in  step  11  along  with  the 
criteria  expressed  by  Eq.  (5.2.31)  to  verify  nonoverlapping  of 
grid  lines. 

(13)  Determine  the  locations  of  the  grid  points  in  the  physical 
domain  at  the  initial  time  level  (n=0)  by  using  Eqs.  (5.2.1)  and 
(5.2.2)  with  the  parameters  evaluated  in  the  previous  steps. 

(14)  Determine  the  coefficients  Ka,  Kb,  and  Kc  by  using  Eq.  (5.2.21) 
with  the  grid  system  obtained  in  step  13.  If  K < Ka,  then  grid 
lines  cluster  near  the  boundaries.  If  Kc  > K > Kb,  then  grid 
lines  stretching  away  from  the  boundaries. 

(15)  Specify  the  maximum  deviation  from  orthogonality  permitted  at 
boundaries  (/S)  . 

(16)  Compare  the  deviation  from  orthogonality  6 determined  from  Eq. 
(5.2.40)  with  0 given  in  step  15.  If  6 < p then  continue  the 
algorithm  from  step  19  otherwise  continue  from  step  17. 

(17)  Determine  the  value  of  coordinate  A rj*  by  using  Eq.  (5.2.41). 

(18)  Determine  the  new  number  of  grid  points  JL*  by  using  Eq. (5.2.44) 
and  repeat  the  algorithm  from  step  3. 

(19)  At  this  point,  use  the  same  procedure  described  by  steps  8 to  14 
in  the  algorithm  for  2-D,  rectangularly  shaped  spatial  domains 
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presented  in  Chapter  IV.  Here,  it  is  noted  that  the  gradient 
used  to  determine  the  weighting  function  in  Eq.  (3.2.5)  for  the 
x'  direction , may  be  different  from  the  one  used  for  the  y- 
direction. 

(20)  Replace  the  time  level  n by  n+1  and  calculate  the  total  time 
elapsed  by  using  Eq.  (3.3.11).  Repeat  steps  11  to  19  until  the 
total  time  elapsed  equals  to  the  duration  of  interest  specified 
in  step  5. 

This  completes  the  description  of  the  solution-adaptive,  grid 
generation  method  for  two-dimensional,  arbitrarily  shaped  spatial 
domains  in  which  the  boundaries  can  move  or  be  stationary. 


CHAPTER  VI 


APPLICATIONS  OF  STRETCHING  FUNCTIONS 

In  this  chapter,  the  technique  developed  in  Chapter  III  is  used 
to  formulate  stretching  functions  for  distributing  grid  points  in 
one -dimensional  spatial  domains.  This  chapter  is  divided  in  three 
sections.  In  the  first  section,  a stretching  function  is  formulated  to 
be  used  with  steady-state  problems  in  which  the  solution  is  known  quali- 
tatively but  not  quantitatively.  This  stretching  function  did  not  have 
to  be  coupled  to  the  solution  because  the  regions  where  grid  points  were 
most  needed  were  known  a priori.  In  the  second  section,  stretching 
functions  were  formulated  for  unsteady  problems  in  which  the  solution  at 
some  specified  time  is  known.  In  these  cases,  it  is  shown  that  the 
stretching  function  should  be  coupled  to  the  solution  to  determine  the 
locations  of  grid  points  as  the  problem  is  being  computed.  Finally,  in 
the  third  section,  the  method  developed  was  used  to  generate  grid  points 
necessary  to  solve  a one -dimensional  initial-value  problem  using  finite- 
difference  methods.  The  solution  obtained  in  this  section  by  using  the 
SAAGG  technique  was  compared  to  the  known  analytical  solution,  as  well 
as  to  the  solution  obtained  by  using  equally  spaced  grid  points. 

iL_l Stretching  Functions  Not  Coupled  to  the  Solution 

The  technique  developed  in  Chapter  III  was  designed  for  formulating 
stretching  functions  which  are  coupled  to  the  solution.  However,  the 
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same  technique  also  can  be  used  for  formulating  stretching  functions 
which  are  not  coupled  to  the  solution.  In  order  to  formulate  such  non- 
coupled  stretching  functions,  it  is  necessary  to  specify  the  weighting 
function  appearing  in  Eq.  (3.1.25)  based  on  the  qualitative  interpre- 
tation of  the  physics  of  the  problem.  An  example  of  a non-coupled 
stretching  function  is  formulated  in  this  section. 

Suppose  that  the  qualitative  interpretation  of  the  problem  being 
analized  suggests  that  grid  points  be  clustered  toward  one  boundary  of 
the  spatial  domain.  As  noted  earlier,  grid  points  will  be  clustered  in 
regions  where  the  weight  function  is  high.  Here,  the  weighting  function 
was  chosen  to  be  the  inverse  of  a parabola  whose  vertex  is  in  the 
vicinity  of  the  boundary  where  the  grid  points  should  be  clustered. 
Figure  6 . 1 shows  this  weighting  function  and  the  boundaries  of  the 
spatial  domain  which  are  at  L^-l  and  L2-L. 

The  mathematical  expression  for  that  weighting  function  is 


C2 

W(x) (6.1.1) 

Ax^  + Bx  + C 

where  the  constant  C2  is  defined  by  Eq.  (3.1.24)  and  the  constants  A,  B 
and  C are  arbitrary  constants.  In  this  problem,  the  x-coordinate  of  the 
vertex  of  the  parabola  appearing  in  the  denominator  of  Eq.  (6.1.1)  was 
set  at  x=l/2 . Thus , the  relationship  between  A and  B can  be  determined 
by  setting  the  derivative  of  the  denominator  in  Eq.  (6.1.1)  at  x-1/2  to 
zero,  yielding 


B - -A 


(6.1.2) 
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Figure  6.1  Weighting  function  for  generating  stretching 
functions  to  cluster  grid  points  toward  one 
boundary  of  the  spatial  domain. 
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By  substituting  Eq. (6.1.2)  into  Eq. (6.1.1),  the  following  weighting 
function  is  obtained: 

C2 

W(x) (6.1.3) 

Ax^  - Ax  + C 


With  the  above  weighting  function,  Eq.  (3.1.25)  gives  the  following 
stretching  function: 


K(x) 


1 + 


- Ax  + C)dx 


x 


(6.1.4) 


Integrating  the  above  equation  gives  the  locations  of  the  grid 
points  after  stretching,  as  follows: 


Ax3  Ax2  A A 

( + Cx  ) - ( + C - 1 ) (6.1.5) 

3 2 3 2 


where  x is  defined  by  Eq.  (3.1.4).  Since  A and  C are  arbitrary 
constants,  the  equation  above  can  be  simplified  by  setting  the  second 
parenthesis  in  Eq.  (6.1.5)  to  zero.  This  gives  the  following 
relationship  between  A and  C: 


A A 

+ C -1  - 0 

3 2 


A 

— (6.1.6) 

6 
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By  substituting  Eq. (6.1.6)  into  Eq. (6.1.5)  and  factoring,  the  following 
equation  is  obtained: 

A 

x - x [ 1 ( x - 0.5  )(  1 - x )]  (6.1.7) 

3 

By  comparing  the  above  equation  with  other  stretching  functions  for 
the  same  purpose  appearing  in  the  literature,  it  was  found  that  Eq. 
(6.1.7)  is  identical  to  the  stretching  function  proposed  by  Vinokur  [24] 
who  developed  a method  for  formulating  stretching  functions  based  on 
truncation  error  analysis.  The  Vinokur's  stretching  function  designed 
to  cluster  grid  points  toward  one  boundary  of  the  spatial  domain  is 
written  below  for  comparison. 

x - x [ 1 - 2(1-B) ( x - 0.5  )(  1 - x )]  (6.1.8) 


where  B is  a constant.  By  comparing  the  above  equation  with  Eq. (6.1.7) 
obtained  by  using  the  technique  developed  in  this  investigation,  it  can 
be  seen  that  the  term  2(1-B)  appearing  in  Eq.  (6.1.8)  is  equivalent  to 
the  term  A/3  appearing  in  Eq.  (6.1.7). 

— Stretching  Functions  Coupled  to  the  Solution 
In  this  section,  three  examples  are  given  to  demonstrate  the  ca- 
pability of  the  method  developed  for  formulating  stretching  functions 
that  are  coupled  to  the  solution  as  the  problem  is  being  computed. 
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2 • 1 — Steep  Gradients  in  One  Side  of  the  Spatial  Domain 

Suppose  that  the  gradient  of  the  solution  for  the  problem  being 
&-n&\yze(l  is  known  at  some  time  level  and  is  given  by 

du 

- 0 . 075 (x  - i) -1.075  (6.2.1) 

dx 


Also,  suppose  that  the  number  of  grid  points  to  be  used  for  solving  this 
problem  is  fixed  at  IL-21  and  that  the  boundaries  of  the  spatial  domain 
are  at  L^-l  and  L2=21. 

By  substituting  Eq.  (6.2.1)  into  Eq.  (3.2.5)  gives  the  weighting 
function  which  is 

or 

W(x)  = [ 1 + | 0 . 075 (x  - l)’1-075!  ] (6.2.2) 


With  the  above  weighting  function,  the  stretching  function  for  this 
problem  is  given  by  Eq.  (3.1.25);  i.e., 


'X 

- SF 

[ 1+ 1 0 . 075 (x  - 1) -1  - 075 | ] dx 

Jl 

1 + 20  

f21 

- SF 

[ 1+ | 0 . 075 (x  - 1) -1-075 | ] dx 

K(x)  = — 

x 


(6.2.3) 


where  SF  is  determined  by  Eq.  (3.2.13). 

The  equation  above  is  solved  numerically  and  the  grid  system  thus 
generated  is  shown  in  Fig.  6.2.  In  that  figure  it  can  be  seen  that  the 


89 


Figure  6.2  Grid  system  generated  for  problems  with  steep 
gradients  in  one  side  of  the  spatial  domain. 

(a)  Grid  system  before  stretching. 

(b)  Grid  system  after  stretching. 
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grid  points  are  clustered  near  the  boundary  where  the  gradient  of  the 
solution  is  high. 

^2-2 — Steep  Gradients  in  the  Middle  of  the  Spatial.  Domain 

Suppose  that  the  gradient  of  the  solution  for  the  problem  being 
analyzed  is  known  at  some  time  level  and  is  given  by  the  following 
equations : 


du 


dx 


0.07(11  - x)-°-93 


1 < x < 11 


(6.2.4) 


du 


dx 


0 . 07 (x  - ll)-°-93 


11  < x < 21 


(6.2.5) 


Also  suppose  that  the  number  of  grid  points  to  be  used  for  solving  this 
problem  is  fixed  at  IL-21  and  the  boundaries  of  the  spatial  domain  are 
at  Lq=l  and  L2=21. 

By  substituting  Eqs . (6.2.4)  and  (6.2.5)  into  Eq.  (3.2.5)  gives  the 
weighting  functions  which  are 


W(x)  - [ 1 + | 0.07(11  - x)'°-93|  ]SF  1 < x < 11 

or 

W(x)  - [ 1 + | 0 . 07 (x  - 11) • 93 | ] 11  < x < 21 


(6.2.6) 


(6.2.7) 


With  the  above  weighting  functions,  the  stretching  function  for  this 
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problem  is  given  by  Eq.  (3.1.25),  i.e,  for  1 < x < 11 


rx 

- SF 

[ 1+ | 0 . 07 (x  - 11) -°- 93 | ] dx 

Jl 

1 + 10  

‘11 

[ 1+ 1 0 . 07 (x  - 11) - 0 • 93 | ]'SFdx 

Jl 

(6.2.8) 

x 


and  for  11  < x < 21 


K(x) 


*x 

[ 1+ 1 0 . 07 (x  - 11) ” 0 * 93 | ]’SFdx 

J 11 

11  + 10  

•21 

[ 1+ | 0 . 07 (x  - 11) “ 0 • 93 | ]'SFdx 

J 11 


x 


(6.2.9) 


where  SF  is  determined  by  Eq.  (3.2.13). 

The  equations  above  are  solved  numerically  and  the  grid  system  thus 
generated  is  shown  in  Fig.  6.3.  In  that  figure  it  can  be  seen  that  the 
grid  points  are  clustered  in  the  middle  of  physical  domain  where  the 
gradient  of  the  solution  is  high. 

2 • ? Steep  Gradients  in  Both  Sides  of  the  Spatial  Domain 

Suppose  that  the  gradient  of  the  solution  for  the  problem  being 
analyzed  is  known  at  some  time  level  and  is  given  by  the  following 
equations : 
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1 2 1 

(b) 


Grid  system  generated  for  problems  with  steep 
gradients  in  the  middle  of  the  spatial  domain. 

(a)  Grid  system  before  stretching. 

(b)  Grid  system  after  stretching. 


Figure 


6.3 
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du 

- 1 • 3 (x  - l)'2-3  1 < x < 11  (6.2.10) 

dx 

du 

- 1.3(21  - x)  "2 • 3 11  < x < 21  (6.2.11) 

dx 

Also  suppose  that  the  number  of  grid  points  to  be  used  for  solving  this 
problem  is  fixed  at  IL-21  and  that  the  boundaries  of  the  spatial  domain 
are  at  L^=l  and  L2~21. 

By  substituting  Eqs . (6.2.10)  and  (6.2.11)  into  Eq.  (3.2.5)  gives 

the  weighting  functions  which  are 


or 

W(x)  - [ 1 + 1 1 . 3 (x  - 1) -2  - 3 | ] 


1 < x < 11 


(6.2.12) 


c p 

W(x)  - [ 1 + | 1 . 3 (21  - x)'2-3|  ] 11  < x < 21 


(6.2.13) 


With  the  above  weighting  functions,  the  stretching  function  for  this 
problem  is  given  by  Eq.  (3.1.25);  i.e.,  for  1 < x < 11 


*x 

[1+| 1. 3 (x  - l)-2-3| j'SFdx 

Jl 

1 + 10  — 

•11 

1+ | 1 . 3 (x  - l)-2-3|]'SFdx 

Jl 


K(x)  - 


x 


(6.2.14) 
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and  for  11  < x < 21 


rx 

- <SF 

[i+ 1 1 . 3 (x  - ir2-3|  dx 

Jl 

11  + 10  

r21 

- SF 

[ 1+ | 1 . 3 (x  - l)-2-3|]  dx 

J 11 

K(x)  - (6.2.15) 

x 

where  SF  is  determined  by  Eq.  (3.2.13). 

The  equations  above  are  solved  numerically  and  the  grid  system  thus 
generated  is  shown  in  Fig.  6.4.  In  that  figure  it  can  be  seen  that  the 
that  the  grid  points  are  clustered  toward  the  two  boundaries  of  the 
spatial  domain  where  the  gradient  of  the  solution  is  high. 

These  three  examples  shown  the  capability  of  the  method  developed  to 
formulate  stretching  functions  so  that  the  grid  points  in  the  spatial 
domain  are  distributed  where  they  are  most  needed. 

6_.,3  Stretching  Functions  for  1-D  Initial -Value  Prohlams 
In  this  section,  the  method  for  formulating  stretching  functions  was 
used  to  generate  a grid  system  for  solving  a one -dimensional  initial- 
value  problem.  The  problem  was  solved  by  using  a finite-difference 
method  and  a grid  system  generated  by  a stretching  function.  The 
solution,  thus  obtained  was  compared  with  the  solution  obtained  by  using 

a grid  system  with  equally  spaced  grid  points,  and  the  exact  known 
solution. 

The  initial-value  problem  to  be  solved  numerically  is  given  by 
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Figure  6.4  Grid  system  generated  for  problems  with  steep 
gradients  in  both  sides  of  the  spatial  domain. 

(a)  Grid  system  before  stretching. 

(b)  Grid  system  after  stretching. 
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au  -1.725 

+ 0.725  x - 0 (6.3.1) 

dx 

subjected  to  the  following  condition: 

u(x-l)  - 1.0  (6.3.2) 

The  number  of  grid  points  used  for  solving  this  problem  is  IL-21  and  the 
boundaries  of  the  spatial  domain  are  L^-l  and  L2-H. 

The  solution  of  this  problem  was  obtained  by  using  the  SAAGG  method 
following  the  steps  of  the  algorithm  presented  in  section  3.3  of  Chapter 
III.  The  first  step  in  that  algorithm  is  to  transform  Eq.  (6.3.2)  so 
that  the  £ is  the  independent  variable  instead  of  x.  The  transformed 
equation  for  this  problem  in  given  by 

du  -1.725 

Cx  + 0.725  [x(0]  - 0 (6.3.3) 

d? 

where  x(£)  is  determined  by  using  Eq . (3.2.2)  with  the  technique  for 

generating  stretching  functions  and  the  metric  coefficient  £x  is  de- 
termined by  using  Eq.  (3.3.8). 

With  the  transformed  equation  given  by  Eq.  (6.3.3),  the  next  step  is 
to  use  a finite-difference  method  to  derive  a finite  difference  equation 
(FDE) . Here,  a first-order-accurate,  backward- difference  approximation 
was  used  and  the  following  FDE  was  obtained  from  Eq.  (6.3.3): 

ui  • ui- 1 c -1.725 

?Xi  + 0.725[x(O]  “ 0 (6.3.4) 

AC 

The  numerical  solution  of  Eqs.  (6.3.1)  and  (6.3.2)  was  obtained  from 
Eq.  (6.3.4)  and  the  results  are  shown  in  Fig.  6.5.  This  figure,  also 
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1 21 

Figure  6.5  Solutions  obtained  for  the  initial -value  problem  given 
by  Eqs . (6.3.1)  and  (6.3.2).  The  unit  for  u is  m/sec  and 
the  unit  for  x is  m. 

(A)  Numerical  solution  obtained  by  using  a grid  system 
with  equally  spaced  grid  points. 

(B)  Numerical  solution  obtained  by  using  stretching 
function  and  the  SAAGG  method. 

(C)  Exact  solution. 
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shows  a comparison  between  the  solution  obtained  by  using  the  SAAGG 
method  developed  in  this  investigation,  the  exact  solution,  and  the 
numerical  solution  obtained  by  using  a grid  system  with  equally  spaced 
grid  points.  In  this  comparison,  it  can  be  seen  that  numerical  errors 
were  reduced  when  the  problem  was  solved  by  using  a stretching  function 
for  redistributing  the  grid  points  in  the  spatial  domain.  Here,  it  is 
noted  that  this  improvement  was  obtained  by  using  the  same  number  of 
grid  points  and  the  same  finite -difference  method.  In  Appendix  C,  the 
computational  efficiency  of  the  SAAGG  method  is  compared  with  the 
computational  efficiency  of  methods  using  equally  spaced  grid  points. 

This  concludes  the  applications  of  the  technique  for  formulating 
stretching  functions  developed  in  Section  3.1  of  Chapter  III.  The 
examples  given  in  this  section  demonstrated  the  feasibility  and  useful- 
ness of  the  technique  investigated.  In  the  next  chapter,  more  examples 
are  given  for  using  this  technique  in  conjunction  with  the  SAAGG  method 
for  solving  one -dimensional  unsteady  problems. 


CHAPTER  VII 


RESULTS  FOR  ONE -DIMENSIONAL  UNSTEADY  PROBLEMS 

In  this  chapter,  the  solution- adaptive , algebraic  grid  generation 
(SAAGG)  method  developed  in  Chapter  III  was  used  with  finite-difference 
methods  to  obtain  numerical  solutions  for  two  one -dimensional  problems. 
The  two  problems  analyzed  were  the  one -dimensional  linearized  inviscid 
Burgers'  equation  and  the  quasi- linear  inviscid  Burgers'  equation.  The 
numerical  solutions  obtained  for  these  two  problems  were  compared  with 
the  exact  solutions. 


7.1  Linearized  Inviscid  Burgers'  Equation 
equation  has  been  used  widely  as  a model  equation  for 
testing  and  comparing  computational  techniques  in  fluid  dynamics.  This 
is  because  it  is  mathematically  simple;  it  describes  some  of  the  most 
important  transport  mechanisms  in  fluid  dynamics;  and  for  certain 
initial  and  boundary  conditions,  it  is  easy  to  obtain  the  exact  solution 
in  analytical  form  [30,31],  The  one-dimensional  Burgers'  equation  is 


du  dF 

+ 

3t  3x 


(7.1.1) 


where  F is  a function  of  u,  and  v is  a constant. 
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The  one -dimensional,  linearized  inviscid  Burgers'  equation  is  ob- 
tained from  Eq.  (7.1.1)  by  neglecting  the  second-order  derivative  term 
and  setting  the  function  F to  cu,  in  which  c is  a constant.  In  this 
section,  it  is  desired  to  use  the  SAAGG  method  along  with  a finite- 
difference  method  to  obtain  numerical  solutions  to  the  one - dimens ional , 
linearized  inviscid  Burgers'  equation  given  below: 

flu  flu 

+ c - 0 (7.1.2) 

3t  Ax 

subjected  to  the  following  initial  and  boundary  conditions: 


u(x,0) 


2(4  - x) 
3 
1 


1 < x < 2.5 


x > 2.5 


(7.1.3) 


u(L]_ , t)  - 2 (7.1.4) 

In  the  equations  above,  c*=lm/sec  and  x varies  from  L^-l  m to  L,2**ll  m. 

In  order  to  use  the  SAAGG  method,  it  is  necessary  to  follow  the 
steps  of  the  algorithm  presented  in  Section  3.3  of  Chapter  III.  The 
first  step  in  the  algorithm  is  to  transform  Eq.  (7.1.2)  so  that  f and  r 
are  the  independent  variable  instead  of  x and  t.  For  the  problem 
described  by  Eqs.  (7.1.2)  to  (7.1.4),  the  transformation  of  the  physical 
domain  in  x-t  coordinate  system  to  the  computational  domain  in  £-r 
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coordinate  system  is  given  by 


£ “ £(x,t)  (7.1.5) 

r - t (7.1.6) 

and  the  partial  derivatives  with  respect  to  x and  t are  related  to  the 
partial  derivatives  with  respect  to  f and  r by 


d 


dx 


d/-  a 

+ 

dx  df 


dr  d 


dx  dr 


d 


dt 


d£  d 

+ 

dt 


dr  d 


dt  dr 


where 


dr 

- 0 and 

dx 


dr 

1 

3t 


(7.1.7) 


(7.1.8) 


(7.1.9) 


Substitution  of  Eqs . (7.1.7)  to  (7.1.9)  into  Eq.  (7.1.2)  gives  the 
following  transformed  equation: 


d£  d d dt  d 

( t ) u + c ( ) u - 0 

dt  d£  dr  dx  d£ 


which  simplifies  to 


102 


3u  3u 

+ dt  + c£x>  “ 0 (7.1.10) 

3r  di 

With  the  governing  equation  transformed,  the  next  step  is  to  use  a 
finite-difference  method  to  derive  a finite-difference  equation  (FDE) . 
Here,  the  Lax-Wendroff  method  [29]  was  used  and  the  following  FDE  was 
derived  from  Eq.  (7.1.10): 

n+1  n Ar  n n 

ui  “ U£  - (£t  + c£x)i  (ui+1  - Ui.x)  + 

2A£ 

, Ar  (2)  (2)  n n n 

< > (ft  + cfx>i  (ui+l  - 2Ui  + ui.!)  (7.1.11) 

2A| 

In  the  equation  above,  the  superscript  (2)  denotes  square  and  the 
superscript  n and  n+1  denote  time  levels.  The  metric  coefficients 
and  are  obtained  from  Eqs . (3.3.8)  and  (3.3.9),  respectively. 

The  next  two  steps  of  the  algorithm  are  to  specify  the  duration  of 
interest  (T)  and  the  number  of  grid  points  (IL)  to  be  employed.  Here,  T 
was  set  equal  to  9 seconds  and  IL  was  set  to  21.  The  remaining  steps  of 
the  algorithm  described  in  Section  3.3  are  straightforward  and  are  not 
repeated  except  for  the  procedure  used  to  evaluate  the  weighting 
function. 

The  weighting  function  appearing  in  Eq.  (3.2.5)  need  to  be  repre- 
sented in  terms  of  the  £-  and  r -coordinates . This  can  be  done  by  using 
Eqs.  (7.1.7)  to  (7.1.9)  as  shown  below: 
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3u  gp 

W(e,r)  - (1  + I fxl)  (7.1.12) 

ae 

By  replacing  3u/3£  in  the  above  equation  by  a second-order-accurate 
m space  central-difference  approximation,  the  finite  difference  form  of 
the  weighting  function  is  obtained,  and  is 

n n 

n ui+l  * ui-l  n qp 

W(|,r)i  = [1  + |( ) CXi  |]  (7.1.13) 

2A£ 

The  weighting  function  in  the  algorithm  described  in  Section  3.3  is 
evaluated  by  using  the  equation  above. 

The  numerical  solution  obtained  from  Eq.  (7.1.2)  to  (7.1.4)  by  using 
Eq.  (7.1.11)  and  the  grid  system  generated  by  the  SAAGG  method  is  shown 
in  Fig.  7.1.  In  that  figure,  it  can  be  seen  that  the  grid  points  in 
physical  domain  are  relocated  after  each  time  step  to  follow  the  steep 
gradient  of  the  problem  which  moves  as  time  progresses. 

The  SAAGG  method  is  more  efficient  than  methods  based  on  grid 
systems  using  equally  spaced  grid  points  because  less  grid  points  are 
needed  to  achieve  the  same  accuracy.  In  this  problem,  for  instance,  the 
number  of  grid  points  would  be  multiplied  by  five  if  a grid  system  using 

equally  spaced  grid  points  was  used  instead  of  the  grid  system  generated 
by  the  SAAGG  method. 

Here,  it  is  noted  that  the  oscillations  near  the  steep  gradient  are 
not  due  to  the  grid  generation  technique  used,  but  are  due  to  the  dis- 
persive nature  of  the  Lax-Wendroff  method  chosen  [23]. 
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Figure  7.1  Solutions  for  the  one - dimens ional , linearized  inviscid 
Burgers'  equation  given  by  Eqs . (7.1.1)  to  (7.1.3).  The 
units  for  x,  u and  t are  m,  m/sec  and  sec,  respectively. 
Numerical  solution  obtained  with  the  grid  system 
generated  by  using  the  SAAGG  method. 

Exact  solution. 
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The  results  given  above  demonstrate  that  the  SAAGG  method  developed 
in  this  investigation  is  able  to  generate  a solution-adaptive  grid 
system  for  solving  one -dimensional  unsteady  problems  in  which  steep 
gradients  can  move  about. 

7.2  Quasi -linear  Inviscid  Burgers'  Equation 
Similar  to  the  linearized,  inviscid  Burgers'  equation,  the  quasi- 
linear  inviscid  Burgers'  equation  is  also  obtained  from  Eq.  (7.1.1)  by 
neglecting  the  second- order  derivative  term  and  setting  the  function  F 
to  u /2.  In  this  section  the  SAAGG  method  is  used  along  with  a finite- 
difference  method  to  obtain  numerical  solutions  for  the  one  - dimens ional , 
quasi-linear  inviscid  Burgers'  equation  given  below: 

2 

3u  3 u 

+ ( ) - 0 (7.2.1) 

at  dx  2 

subjected  to  the  following  initial  and  boundary  conditions: 


u(x , 0) 


2 

1 


sin[ (x  - 2)] 

2 


1 < x < 3 
x > 3 


(7.2.2) 


ux(Li,t)  - 0 


(7.2.3) 


where  x varies  between  L^-l  m and  L2-8 . 5 m. 
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In  order  to  use  the  SAAGG  method,  it  is  necessary  to  follow  the 
steps  of  the  algorithm  presented  in  Chapter  III  and  already  used  in 
Section  7.1.  The  first  step  in  the  algorithm  is  to  transform  Eq. (7.2.1) 
so  that  £ and  r are  the  independent  variable  instead  of  x and  t.  For 
the  problem  describe  by  Eqs.  (7.2.1)  to  (7.2.3),  the  transformation  of 
the  physical  domain  in  x-t  coordinate  system  to  the  computational 
domain  in  £-r  coordinate  system  is  given  by  Eqs.  (7.1.5)  to  (7.1.6),  and 
the  partial  derivatives  with  respect  to  x and  y are  related  to  the 
partial  derivatives  with  respect  to  £ and  r by  Eqs.  (7.1.7)  to  (7.1.9). 
Thus,  substitution  of  Eqs. (7.1.5)  to  (7.1.9)  into  Eq.  (7.2.1)  gives  the 
following  transformed  equation. 

2 

3u  3u  3 u 

+ ft + fx < ) - 0 (7.2.4) 

dr  3£  3£  2 


With  the  governing  equation  transformed,  the  next  step  is  to  use  a 
finite-difference  method  to  derive  a finite-difference  equation  (FDE) . 
Here,  the  Lax-Wendroff  method  [29]  was  used  and  the  following  FDE  was 
derived  from  Eq.  (7.2.4): 


n+1  n 

ui  - Ui  - 


At 

ft (ui+l 

2A£ 


n Ar  (2)  (2)  n 

ui-l)  ' fx (ui+l  * uj,.!)  + 

4Ax 


(2) 

fx  Ar  (2)  n (2) 

( ) [(Ui  + Ui+1)  (ui+1  - 

8 A£ 


(2)  n n (2) 

Ui  ) - (Ui  + Ui_i)  (Ui  - 


(7.2.5) 


(2)  n 
ui-l>  ] 
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In  the  equation  above,  the  superscript  (2)  denotes  square  and  the 
superscript  n and  n+1  denote  time  levels.  The  metric  coefficients  £x 
and  are  obtained  from  Eqs . (3.3.8)  and  (3.3.9),  respectively. 

The  next  two  steps  of  the  algorithm  are  to  specify  the  duration  of 
interest  (T)  and  the  number  of  grid  points  (IL)  to  be  employed.  Here,  T 
was  set  equal  to  6 seconds  and  IL  was  set  to  31.  The  same  procedure 
used  in  Section  7.1  to  evaluate  the  weighting  function  in  terms  of  the 
£*  and  r-coordinates  is  used  in  this  problem.  The  remaining  steps  of 
the  algorithm  described  in  Chapter  III  are  straightforward  and  are  not 
repeated  here. 

The  numerical  solution  obtained  from  Eq.  (7.2.1)  to  (7.2.3)  by  using 
Eq.  (7.2.5)  and  the  grid  system  generated  by  the  SAAGG  method  are  shown 
in  Fig.  7.2.  In  this  figure,  it  can  be  seen  that  the  grid  points  in 
physical  domain  are  relocated  after  each  time  step  to  follow  the  steep 
gradient  of  the  problem  which  moves  as  time  progresses.  This  problem 
differs  from  the  problem  solved  in  Section  7.1  because,  in  this  case, 
the  maximum  gradient  in  the  solution  becomes  steeper  at  each  time  step. 
Here,  it  is  noted  that  the  same  remarks  on  the  advantages  of  the  SAAGG 
method,  made  when  solving  the  linearized  inviscid  Burgers'  equation,  are 
valid  for  this  problem. 

To  further  demonstrate  the  usefulness  of  the  SAAGG  method,  it  was 
used  along  with  the  Lax-Wendroff  method  to  obtain  solutions  for  the  one- 
dimensional,  quasi-linear  inviscid  Burgers'  equation  defined  by  Eq. 
(7.2.1),  subjected  to  another  initial  condition.  This  initial  condition 
is  given  by  the  following  equation: 
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Figure  7.2  Solution  of  the  one - dimens ional , quasi-linear  inviscid 
Burgers'  equation  given  by  Eqs.  (7.2.1)  and  (7.2.3) 
obtained  with  the  grid  system  generated  by  using  the 
SAAGG  method.  The  units  for  x,  u and  t are  m,  m/sec 
and  sec,  respectively. 
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Figure  7.3  Solution  of  the  one - dimens ional , quasi-linear  inviscid 
Burgers'  equation  given  by  Eqs.  (7.2.1)  and  (7.2.6) 
obtained  with  the  grid  system  generated  by  using  the 
SAAGG  method.  The  units  for  x,  u and  t are  m,  m/sec 
and  sec,  respectively. 
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It 

u(x,  0)  - 1.25  + 0 . 75sin[ (x  + 0.5)]  x 2:  1 (7.2.6) 

3 

where  x varies  from  L^— 1 m and  L2~13  m,  the  number  of  grid  points  is 
IL-49  and  the  other  conditions  and  specifications  remain  the  same. 

By  using  the  initial  condition  given  above,  the  capability  of  the 
SAAGG  method  for  generating  grid  systems  to  follow  steep  gradients  in 
more  than  one  region  in  the  physical  domain  is  demonstrated. 

The  numerical  solution  obtained  from  Eq.  (7.2.1)  subjected  to  the 
initial  condition  above  and  the  grid  generated  by  the  SAAGG  method  are 
shown  in  Fig.  7.3.  In  this  figure,  it  can  be  seen  that  the  grid  points 

in  physical  domain  are  relocated  after  each  time  step  to  follow  the 
steep  gradients  in  the  physical  domain.  These  results  demonstrate  the 
capability  of  the  SAAGG  method  to  generate  grid  systems  to  follow  more 
than  one  steep  gradient. 

The  examples  above  conclude  the  demonstrations  of  the  SAAGG  method 
developed  in  this  investigation  for  solving  one -dimensional  unsteady 
problems . 


CHAPTER  VIII 


RESULTS  FOR  TWO-DIMENSIONAL  UNSTEADY  PROBLEMS 

In  this  chapter,  the  solution-adaptive,  algebraic  grid  generation 
(SAAGG)  method  developed  in  Chapter  IV  and  V was  used  with  finite- 
dif ference  methods  to  obtain  numerical  solutions  for  two-dimensional 
problems.  This  chapter  is  divided  in  two  sections.  In  the  first 
section,  the  SAAGG  method  was  used  to  solve  two  two-dimensional  problems 
with  rectangularly  shaped  spatial  domains.  In  the  second  section,  the 
SAAGG  method  was  used  in  conjunction  with  the  Two -Boundary  Technique  to 
generate  a grid  system  in  a convergent-divergent  spatial  domain. 

iLl Problems  With  Rectansularlv  Shaped.  Spatial  Domain 

In  this  section,  two  problems  with  rectangularly  shaped  spatial 
domain  were  solved  by  using  the  SAAGG  method  with  finite-difference 
methods.  The  two  problems  analyzed  were  the  two-dimensional,  unsteady 
problems  of  plane  wave  propagation  and  the  propagation  of  concentric, 
circular  waves  in  the  radial  direction. 

8.1.1  Plane  Wave  Propagation 

The  two-dimensional,  quasi-linear  inviscid  Burgers'  equation  was 
used  with  appropriate  initial  and  boundary  conditions  to  model  a plane 
wave  propagating  in  a direction  45  degrees  from  the  x-axis.  These 
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equations  are  given  below: 


2 2 
3u  3 u 3 u 

+ ( )+  ( ) - 0 (8.1.1) 

at  ax  2 ay  2 


tt(x  + y) 

u(x,y , 0)  - 2 + sin  (8.1.2) 

4 

ux(Ll . , t)  - 0 (8.1.3) 

ux  ( ^1  > ^*1 1 b ) — 0 (8.1.4) 

where  x varies  between  L^-l  m and  L,2“7  m,  and  y varies  between  Lj-l  m 
> 

and  L-2-7  m. 

To  use  the  SAAGG  method,  it  is  necessary  to  follow  the  steps  of  the 
algorithm  presented  in  Chapter  IV.  The  first  step  is  to  transform  Eq. 
(8.1.1)  so  that  £,  r/  and  r are  the  independent  variable  instead  of  x,  y 
and  t.  For  the  problem  described  by  Eqs . (8.1.1)  to  (8.1.4),  the 

transformation  of  the  physical  domain  in  x-y-t  coordinate  system  to  the 
computational  domain  in  £-r;-r  coordinate  system  is  given  by  Eq.  (2.3) 
and  the  partial  derivatives  with  respect  to  x,  y and  t are  related  to 
the  partial  derivatives  with  respect  to  (,  i|  and  r by  Eq.  (2.4).  Thus, 
substitution  of  Eqs.  (2.3)  and  (2.4)  into  Eq.  (8.1.1)  gives  the  fol- 
lowing transformed  equation: 
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3u 

a? 


3u  3u 

+ r,t + 

dr,  dr 


2 

d u 

+ <£x  + Cy) — ( — ) 

3£  2 


+ 


2 

d u 

(*7x  + Vy) ( ) - o (8.1.5) 

dr,  2 


The  equation  above  can  be  written  more  compactly  as  follows: 


3u  3F  3G 


dr  d£  dr. 


(8.1.6) 


where 


F 


2 

u 


£tu  + (£x 


+ Cv)- 


(8.1.7) 


2 

u 

G - r?tu  + (»7X  + fjy) 

2 


(8.1.8) 


With  the  governing  equation  transformed,  the  next  step  in  the 
algorithm  is  to  use  a finite  difference  method  to  derive  a finite- 
difference  equation  (FDE) . Here,  the  Lax-Wendroff  method  [29]  was  used 
and  the  following  FDE  was  derived  from  Eq . (8.1.6): 
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n n n n 

n+1  n Ar  Fi+1j  - F^ij  A r Glfj+1  - Gitj_  ! 

ui,  j - ui  . j 

A£  2 At?  2 


1 At  (2)  n n n n 

tAi+l/2(Fi+l,j  * Fi , j ) * Ai-l/2(Fi,j  - Fi-l,j)]  + 


1 , Ar  (2)  n n n n 

( > [Bi+l/2(Gi,j+l  ' Gi,j)  - Bj - 1/2 <Gi , j ’ Gi,j-1>] 

2 Ar/ 


where 


n 

Fi, 


n 

ui 


+ (£x  + £v) 


n(2> 

ui.j 


n 


n 


;i,j  “ »Jt  ui , j + <»7x  + »?y)' 


n(2) 

ui.j 


n n 

Ai±l/2  “ £t  + “ + £y)  (ui , j + «i±l,j) 


n n 

Bj±l/2  ' Vt  + — C»?x  + »?y)(ui,j  + ui,j±l) 


(8.1.9) 


(8.1.10) 


(8.1.11) 


(8.1.12) 


(8.1.13) 
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n 

Fi±l,j 


n 

“ £t  ui±l , j + (Cx 


+ Cv) 


n(2> 

ui±l,j 


(8.1.14) 


Gi,j±l 


n 


“ »?t  ui, j±l  + (»?x  + Vy) 


n(2) 

ui,j±l 


(8.1.15) 


In  the  equation  above,  the  superscript  (2)  denotes  square  and  the 
superscripts  n and  n+1  denote  time  levels.  The  metric  coefficients 

fy>  *?x > *7y>  £t*  *?t>  anc^  Tt  are  obtained  from  Eqs.  (2.6)  to 
(2.12),  respectively. 

The  next  two  steps  of  the  algorithm  are  to  specify  the  duration  of 
interest  (T)  and  the  number  of  grid  points  (IL)  and  (JL)  to  be  employed. 
Here,  T was  set  equal  to  6 seconds  and  IL  and  JL  were  both  set  to  13. 
The  remaining  steps  of  the  algorithm  described  in  Section  4.1  are 
straightf oxrward  and  are  not  repeated  here . The  weighting  function  in 
the  f -direction  is  obtained  from  Eq.  (7.1.13).  The  weighting  function 
in  the  ^-direction  is  also  obtained  from  Eq.  (7.1.13)  by  replacing  the 
coordinate  f by  r/ . 

The  numerical  solution  obtained  from  Eqs.  (8.1.1)  to  (8.1.4)  by 
using  Eq.  (8.1.9)  and  the  grid  system  generated  by  the  SAAGG  method  are 
shown  in  Figs.  8.1  to  8.5,  where  the  unit  for  x and  y is  m and  the  unit 
for  u depends  on  the  variable  being  analyzed.  In  these  figures,  it  can 
be  seen  that  the  grid  points  in  the  two-dimensional  physical  domain  are 
relocated  after  each  time  step  to  follow  the  steep  gradient  of  the  plane 


wave  as  it  travels. 
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Figure  8.1  Solution  for  the  propagation  of  a plane  wave  obtained 
by  using  the  grid  system  generated  by  the  SAAGG  method 
at  t-1  sec. 
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Figure  8.2  Solution  for  the  propagation  of  a plane  wave  obtained 

by  using  the  grid  system  generated  by  the  SAAGG  method 
at  t-1.5  sec. 


118 


SIDE  VIEW 


7 

SAAGG  2-D  GRID 


Figure  8 . 3 Solution  for  the  propagation  of  a plane  wave  obtained 
by  using  the  grid  system  generated  by  the  SAAGG  method 
at  t-2  sec. 
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Figure  8.4  Solution  for  the  propagation  of  a plane  wave  obtained 
by  using  the  grid  system  generated  by  the  SAAGG  method 
at  t-3  sec. 
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Figure  8.5  Solution  for  the  propagation  of  a plane  wave  obtained 
by  using  the  grid  system  generated  by  the  SAAGG  method 
at  t-6  sec. 
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The  remarks  made  in  Section  7.1  for  the  one -dimensional  problems  are 
also  valid  for  the  two-dimensional  problems  discussed  in  this  section. 
Again,  using  the  grid  system  generated  by  the  SAAGG  method  is  more 
efficient  than  using  a grid  system  with  equally  spaced  grid  points  to 
achieve  the  same  desired  accuracy  in  the  solution.  In  this  problem,  for 
instance,  the  number  of  grid  points  would  be  multiplied  by  sixteen  if  a 
grid  system  with  equally  spaced  grid  points  was  used  instead  of  the  grid 
system  generated  by  the  SAAGG  method.  In  Appendix  C,  the  SAAGG  method 
is  compared  with  an  equally  spaced  method  in  terms  of  efficiency. 

As  was  noted  for  one -dimensional  problems,  the  oscillations  along 
the  steep  gradient  are  not  due  to  the  grid  generation  technique  used, 
but  to  the  dispersive  nature  of  the  Lax-Wendroff  method  chosen  [23]. 
Here,  it  is  noted  that  the  steep  gradient  tends  to  smear  near  the 
boundaries  after  several  time  steps.  This  is  due  to  the  first-order- 
accurate  finite-difference  approximation  used  to  compute  the  solution  at 
the  boundaries . 

The  results  obtained  in  this  section  demonstrate  that  the  SAAGG 
method  can  generate  solution-adaptive  grid  systems  for  analyzing  two- 
dimensional  unsteady  problems  with  steep  gradients  which  move  about  as 
time  progresses. 

8.1.2  Circular  Wave  Propagation 

The  two-dimensional,  quasi-linear  inviscid  Burgers'  equation  was 
used  with  appropriated  initial  and  boundary  conditions  to  model  the 
propagation  of  concentric  circular  waves  in  the  radial  direction.  These 
equations  are  given  below: 
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2 2 
flu  flu  3 Eu 

+ ( ) + ( ) “ 0 (8.1.16) 

3t  3x  2 fly  2 


9 9 1/2 

tt[2(x2  + y2)] 

u(x ,y , 0)  - 2 + sin  (8.1.17) 

4 


ux(Li,Li,t)  - 0 (8.1.18) 


ux(Ll>Li,t)  - 0 (8.1.19) 


where 

y - L[ 

E - (8.1.20) 

x - L]^ 

In  the  above  equations,  x varies  between  L^-l  m and  L2~7  m,  and  y 
varies  between  l{-1  m and  L2~7  m. 

To  use  the  SAAGG  method,  it  is  necessary  to  follow  the  steps  of  the 
algorithm  presented  in  Chapter  IV.  The  first  step  in  the  algorithm  is 
to  transform  Eq.  (8.1.16)  so  that  £ , r\  and  r are  the  independent 
variable  instead  of  x,  y and  t.  For  the  problem  describe  by  Eqs . 
(8.1.16)  to  (8.1.20),  the  transformation  of  the  physical  domain  in  x-y-t 
coordinate  system  to  the  computational  domain  in  coordinate  system 

is  given  by  Eq.  (2.3)  and  the  partial  derivatives  with  respect  to  x,  y 
and  t are  related  to  the  partial  derivatives  with  respect  to  £ , r)  and  r 
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by  Eq.  (2.4).  Thus,  the  substitution  of  Eqs.  (2.3)  and  (2.4)  into  Eq. 
(8.1.16)  gives  the  following  transformed  equation: 


3u  3u  3u  d u 

ft + Vt + + (£x  + Efy) ( ) + 

df  dr,  dr  d£  2 


2 2 

d u u 

(»?x  + Erly) ( ) + 0 (8.1.21) 

dr,  2 2 (x  - L].) 


The  equation  above  can  be  written  more  compactly  as  follows: 


3u  dF  3G 

S (8.1.22) 

dr  dr. 


where 


2 

u 

F - ft«  + (fx  + E^y) 

2 


(8.1.23) 


2 

u 

»?tu  + (*Jx  + Ef?y) 


S 


2 

u 


2(x  - Li) 


(8.1.24) 


(8.1.25) 


With  the  governing  equation  transformed,  the  next  step  in  the 
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algorithm  is  to  use  a finite  difference  method  to  derive 
difference  equation  (FDE) . Here,  the  Lax-Wendroff  method  [29] 
and  the  following  FDE  was  derived  from  Eq.  (8.1.22): 


n+1  n 

ui,j  “ ui,j 


n n n n 

Ar  F1+1,j  ~ Fj.xj  Ar  G1(j+1  - Gjj.! 

2 At)  2 


1 t Af  , n n n n 

2 ( A?  tAi+l/2 (Fi+1 , j - Fi,j)  - Ai-1/2(Fii j - Fi-l,j)j  + 


1 . Ar  (2)  n n n 

( ) [Bi+l/2(Gi,j+l  * Gi,j)  * Bj-l/2(Gi,j 

l Arj 


G?j-i)]  + 


Aru 


2(x  - LX) 


where 


n 


n 


ri,j  “ Ct  ui , j + (£x  + Ei£y) 


n(2> 

ui.j 


(2) 


n n 

3i,j  “ * t ui,j  + (»7x  + Ei»?y) 


n 

ui.j 


f inite- 
was  used 


(8.2.26) 


(8.1.27) 


2 


(8.1.28) 
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n 


n 


Ai±l/2  “ £t  + — (Cx  + Ei£y) (ui , j + ui±l,j) 


(8.1.29) 


1 

Bj±l/2  - »?t  + ~ 


(»?x  + Ei^y)  (ui , j + ui,j±l) 


(8.1.30) 


n 

Fi±l,j 


n 


” ^t  Ui±l, j + (^X  + Ei^y) 


n(2) 

ui±l,j 


(8.1.31) 


n n ui , j±l 

Gi,j±l  “ »»t  ui, j±l  + (»?x  + Ei^y) (8.1.32) 

2 

In  the  equation  above,  the  superscript  (2)  denotes  square  and  the 
superscripts  n and  n+1  denote  time  levels.  The  metric  coefficients  £x, 
£y.  *7x  > *7y>  £t>  *?t>  and  Tt  are  obtained  from  Eqs . (2.6)  to  (2.12). 

The  next  two  steps  of  the  algorithm  are  to  specify  the  duration  of 
interest  (T)  and  the  number  of  grid  points  (IL)  and  (JL)  to  be  employed. 
Here,  T was  set  equal  to  6 seconds  and  IL  and  JL  were  set  equal  to  13. 
The  remaining  steps  of  the  algorithm  described  in  Section  4.1  are 
straightforward  and  are  not  repeated  here.  The  weighting  function  in 
the  ^-direction  is  obtained  from  Eq.  (7.1.13).  The  weighting  function 
tt*  the  ^-direction  is  also  obtained  from  Eq.  (7.1.13)  by  replacing  the 
coordinate  ( by  i;. 

The  numerical  solution  obtained  from  Eq.  (8.1.16)  to  (8.1.20)  by 
using  Eq.  (8.1.26)  and  the  grid  system  generated  by  the  SAAGG  method  are 
shown  in  Figs.  8.6  to  8.8.  In  these  figures  it  can  be  seen  that  the 
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PERSPECTIVE  VIEWS 

Figure  8.6  Solution  for  the  propagation  of  concentric  circular 
waves  obtained  by  using  the  grid  system  generated  by 
the  SAAGG  method  at  t=2  sec . 
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Figure  8.7  Solution  for  the  propagation  of  concentric  circular 
waves  obtained  by  using  the  grid  system  generated  by 
the  SAAGG  method  at  t-4  sec. 
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Figure  8.8  Solution  for  the  propagation  of  concentric  circular 
waves  obtained  by  using  the  grid  system  generated  by 
the  SAAGG  method  at  t-6  sec . 
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grid  points  in  the  two-dimensional  physical  domain  are  relocated  after 
each  time  step  by  the  SAAGG  method  to  follow  the  steep  gradient  of  the 
circular  wave  travelling  in  the  radial  direction.  This  example  conclu- 
des the  demonstration  of  the  SAAGG  method  developed  in  Chapter  IV  for 
problems  with  rectangularly  shaped  spatial  domain. 

iL-2 Problems  With  Arbitrarily  Shaped.  Spatial  Domains 

In  this  section,  the  solution-adaptive,  algebraic  grid  generation 
(SAAGG)  method  developed  in  Chapter  V is  used  with  the  Two-Boundary 
Technique  to  generate  a grid  system  in  a convergent -divergent  spatial 
domain. 

In  order  to  use  the  SAAGG  method,  it  is  necessary  to  follow  the 
steps  of  the  algorithm  presented  in  Section  5.3  of  Chapter  V.  Here,  the 
steps  1,  2 and  4 of  the  algorithm  are  not  applied  because  it  is  desired 
to  generate  a grid  system  in  a spatial  domain  for  a problem  in  which  the 
solution  is  known  at  a given  time  level.  In  step  3,  the  number  of  grid 
lines  were  specified  as  IL=21  and  JL=13,  and  in  step  5,  the  correspon- 
dence between  the  boundaries  of  the  spatial  domain  and  the  computational 
domain  is  shown  in  Fig.  8.9. 

In  the  next  step  of  the  algorithm,  the  coordinates  of  the  two 
boundaries  AD  and  BC  shown  in  Fig.  8.9  and  expressed  by  Eqs . (5.1.1.) 
to  (5.1.4)  are  given  by 

xl(!.0)  - (L2  - L1)|  + L;l  (8.2.1) 


yi(£,0)  = V - (V  - 1)(2£  - l)2 


(8.2.2) 
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(b) 


Figure  8.9  Correspondence  between  the  boundaries  of  the  physical 
domain  and  the  computational  domain. 

(a)  convergent -divergent  shaped  physical  domain  (x-y-t). 

(b)  computational  domain 


131 


X2(C,1)  “ (L2  - L]_)|  + Li  (8.2.3) 

y2(e,D  - (V  + Dm)  - (V  - 1)(2£  - l)2  (8.2.4) 

where  L]_  - 1 cm  and  L2  « 9 cm,  V-2.25  cm  and  Dm-2.5  cm. 

The  partial  derivatives  appearing  in  Eqs . (5.2.1)  and  (5.2.2)  are 

obtained  from  Eqs.  (8.2.1)  to  (8.2.4)  and  are 

A2  - Ai  = L2  - Li  - 8 (8.2.5) 

B2  — BX  = 4(V  - 1)(2£  - 1)  - 5 (2£  - 1)  (8.2.6) 

With  the  partial  derivatives  obtained  and  the  coordinates  Q and  P 
obtained  from  Eqs.  (5.2.7)  to  (5.2.8),  the  coefficients  are  deter- 
mined from  Eqs.  (5.2.11)  and  (5.2.12)  for  each  grid  line  and  are 
given  by 

Q - xi 

Ki  - i-2,3 IL-1  (8.2.7) 

(V-l)(2fi-l) 

These  coefficients  are  compared  with  the  criteria  for  nonoverlapping 
evaluated  by  using  Eq.  (5.2.31)  with  D-Dm=2 . 5 cm  and  A-L2-L]-8  cm;  i.e., 

0 < Ki  < 0.937  (8.2.8) 


At  this  point,  the  initial  grid  system  can  be  determined  by  using 
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Eqs.  (5.2.1)  and  (5.2.2)  with  the  coefficients  given  by  Eq.  (8.2.5), 
the  partial  derivatives  , A2 , and  B2  given  by  Eqs.  (8.2.5)  and 
(8.2.6),  and  the  coefficients  h]_,  h2 , I13  and  114  given  by  Eqs.  (5.1.7)  to 
(5.1.10).  With  this  initial  grid  system  the  coefficients  Ka,  Kb  and  Kc 
can  be  determined  by  using  Eq.  (5.2.21)  and  are 


Ka 

« 0. 

.316 

(8.2.9) 

Kb 

« 0. 

,567 

(8.2 

.10) 

Kc 

* 0. 

,937 

(8.2. 

.11) 

The  coefficients  determined  above  are  compared  with  the  coefficients  Kj[ 
obtained  from  Eq.  (8.2.7)  to  ensure  that  every  grid  line  is  clustering 
toward  the  boundaries  by  using  the  following  inequality: 

0.567  < Ki  < 0.937  (8.2.12) 

In  the  next  step  of  the  algorithm,  the  deviation  from  orthogonality 
is  specified  as  0—0.9  degrees.  This  value  is  compared  to  the  values  of 
0 i obtained  by  using  Eq.  (5.2.40).  In  this  problem,  every  6 was  found 
to  be  less  than  0.  Thus,  the  algorithm  proceeds  from  step  19. 

In  step  19  of  the  algorithm,  the  gradient  of  velocity  Uy  and  the 
gradient  of  pressure  Px  were  used  with  Eq.  (3.2.5)  to  determine  the 
weighting  functions  W(y)  and  W(x)  , respectively,  to  be  used  with 
the  SAAGG  method  for  determining  the  coordinates  of  the  grid  points  in 
the  spatial  domain. 

The  grid  system  generated  by  using  the  SAAGG  method  is  shown  in  Fig. 
8.10  for  two  different  profiles  of  pressure  inside  the  spatial  domain. 
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Figure  8.10  Grid  system  generated  by  using  the  SAAGG  method 
inside  a convergent -divergent -shaped  spatial 
domain  for  two  different  profiles  of  pressure. 

(a)  Mach  number  1.5. 

(b)  Mach  number  2.0. 
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In  this  figure,  it  can  be  seen  that  the  grid  points  in  the  spatial 
domain  are  clustered  to  regions  of  steep  gradients;  i.e.,  near  the 
region  of  discontinuity  of  pressure  in  the  x-direction  and  near  the 
region  of  high  gradients  of  velocity  at  boundaries  in  the  y-direction. 
Here,  it  is  noted  that  the  SAAGG  method  was  able  to  determine  the 
coordinates  of  grid  points  independently  in  each  coordinate  direction  by 
using  different  weighting  functions,  namely,  the  pressure  gradient  in 
the  x-direction  and  the  velocity  gradient  in  the  y-direction. 

The  examples  above  conclude  the  demonstration  of  the  SAAGG  method 


developed  in  this  investigation  for  generating  solution-adaptive  grid 
systems  in  arbitrarily  shaped  spatial  domains. 


CHAPTER  IX 


CONCLUSIONS  AND  RECOMMENDATIONS 

In  this  investigation,  a solution- adaptive , algebraic  grid  genera- 
tion (SAAGG)  method  was  developed  which  can  generate  non-uniform  grid 
systems  needed  by  finite-difference  algorithms  to  obtain  numerical 
solutions  for  fluid  flow  problems . The  grid  systems  generated  by  the 
SAAGG  method  are  designed  especially  for  unsteady  problems  with 
disparate  length  scales  in  different  parts  of  the  spatial  domain  which 
can  move  about  in  an  unpredictable  manner. 

The  SAAGG  method  is  based  on  a new  class  of  stretching  functions 
developed  in  this  investigation  which  are  derived  from  a variational 
principle.  These  stretching  functions  couple  the  distribution  of  the 
grid  points  in  the  physical  domain  to  the  numerical  solution  being 
computed.  That  technique  for  formulating  stretching  functions  also  can 
be  used  for  formulating  stretching  functions  which  are  not  coupled  to 
the  solution  for  use  with  any  non-adaptive  grid  generation  method.  When 
using  stretching  functions  which  are  coupled  to  the  solution  for 
designing  the  SAAGG  method,  parameters  controlling  the  smoothness  and 
orthogonality  of  the  grid  system  were  added  to  the  grid  generation 
process . This  makes  the  SAAGG  method  more  useful  in  generating  grid 
systems  to  be  used  with  finite-difference  methods. 
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The  SAAGG  method  was  applied  to  several  problems  and  these  applica- 
tions demonstrated  that  the  SAAGG  method  can  redistribute  the  grid 
points  in  the  physical  domain  where  they  are  most  needed  to  reduce 
numerical  errors.  By  placing  the  grid  points  where  most  needed  the 
total  number  of  grid  points  necessary  to  obtain  solutions  of  a desired 
accuracy  can  be  reduced  considerably  and  this  enhance  computational 
efficiency  as  shown  in  Appendix  C. 

As  with  any  other  method,  improvements  and  extensions  can  be  made. 
For  the  SAAGG  method  developed  is  this  investigation,  the  following 
improvement,  extension  and  aplication  are  suggested  in  future  works: 

1.  An  improvement  in  the  technique  developed  for  formulating 
stretching  functions  so  that  smoothness,  orthogonality  and 
clustering  can  all  be  accounted  for  by  a variational  principle. 

2.  An  extension  of  the  SAAGG  method  so  that  it  can  generate 
solution- adaptive  grid  systems  in  3-D  spatial  domains. 

3 . An  application  of  the  SAAGG  method  so  that  it  can  be  combined 
with  the  Four-Boundary  Technique  to  generate  solution -adaptive 
grid  systems  in  which  four  boundaries  of  the  spatial  domain  can 
be  account  in  the  formulation. 


APPENDIX  A 


DERIVATION  OF  THE  EULER- LAGRANGE  EQUATION 

The  Euler -Lagrange  equation  is  a partial  derivative  equation  which 
once  solved  leads  to  a function  y(x)  which  is  a solution  for  an  integral 
problem  of  minimization.  The  problem  is  to  determine  a function  y(x)  so 
that  an  integral  between  two  limits  xq  and  X]^  is  a minimum.  This 
integral  is 


I 


■*1 

F(x,y,y')  dx 

Jx0 


(A.l) 


where  F is  a functional  formed  by  the  independent  variable  x,  the 
dependent  variable  y and  its  first  derivative  y' . 

By  assuming  that  y(x)  is  the  solution  for  the  problem,  the  fol- 
lowing substitution  can  be  made  in  Eq.  (A.l): 


y(x)  - y(x)  + a g(x)  (A. 2) 

where  a is  a constant  between  (— < »,+»)  and  g(x)  is  a function  so  that 
g(xO)“g(xl)“0 • Then,  Eq.  (A.l)  can  be  written  as  follows: 


I(«) 


■*1 

F(x,y  + “g,y'+  ag')  dx 

Jx0 


(A. 3) 
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When  a - 0 the  integral  above  is  a minimum  because  in  this  case 
y(x)  - y(x)  is  the  solution  for  the  problem.  Then  the  variation  of  Eq. 
(A. 3)  is  zero ; i . e . , 


3 I 

51(a)  = - 0 (A. 4) 

3a 


The  substitution  of  Eq.  (A. 3)  into  Eq.  (A. 4)  yields 


51(a) 


•*1  a 

[F(x,y ,y ' ) ] 6a  dx 

. xq  dcc 


(A. 5) 


In  the  above  equation  the  partial  derivative  can  be  evaluate  as  follows: 


5F  3F  3y  3F  3y'  3F  3F 

“ I + g(x)  + g'(x)  (A. 6) 

3a  3y  a a ay'  doc  g y 9y> 


Substitution  of  Eq.  (A. 6)  into  Eq.  (A. 5)  yields 


51(a) 


**1  3F  3 F 

-[ g(x)  + g'(x)]  6a  dx  - 0 

J xq  3x  3y ' 


(A. 7) 


In  Eq.  (A. 7)  the  second  term  can  be  integrated  by  parts,  as  follows: 


^1  3F 

g' (x)dx 

. xo  ay' 


3F  xi 

[ g(x)] 

8y'  x0 


d 3F 

g(x) ( ) 

dx  3y ' 


dx 


x0 


(A. 8) 
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The  first  term  on  the  right  hand  side  of  Eq.  (A. 8)  is  zero  because 
g(*0)=g(x].)«0  by  hypothesis.  Then,  the  substitution  of  Eq.  (A.  8)  into 
Eq.  (A. 7)  yields 


SI(«) 


**1  3F  d 3F 

g(x)  [ - ( )]  6<x  dx  - 0 

Jxq  dy  dx  dy' 


(A. 9) 


The  above  equation  is  satisfied  when  the  term  inside  the  brackets  is 
zero ; i . e . , 


3F  d 3F 

( ) - 0 (A.  10) 

3y  dx  dy' 


The  above  equation  is  known  as  the  Euler -Lagrange  equation  related 


to  the  integral  minimization  problem  of  Eq.  (A.l). 


APPENDIX  B 


EXACT  SOLUTION  OF  THE  LINEARIZED  INVISCID  BURGERS'  EQUATION 


The  linearized  inviscid  Burgers'  equations  was  presented  in  Section 
7 . 1 and  it  is 


3u  3u 

+ c - 0 (B.l) 

5t  dx 


where  c is  a positive  constant. 

The  exact  solution  to  Eq.  (B.l)  can  be  obtained  by  assuming 
continuous  solutions  u-u(t,x),  so  that  the  total  derivative  of  u to  the 
variable  t exists;  i.e., 


du  3u  du  dx 

_ + 

dt  dt  dx  dt 

By  letting 

dx 

- c 

dt 


the  Eq.  (B.2)  can  be  written  as  follows: 


du  du  du 

- + c 

dt  dt  dx 


(B.2) 


(B . 3) 


(B.4) 
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The  substitution  of  Eq.  (B.4)  into  Eq.  (B.l)  gives  the  following 
ordinary  differential  equation  (ODE) . 

du 

0 (B . 5) 

dt 

The  above  equation  implies  that  along  any  characteristic  curve 
defined  by  Eq.  (B.3),  Eq.  (B.l)  becomes  an  ODE. 

The  integration  of  Eq.  (B.3)  gives 


*x 

dx 


X - | + Ct  (B . 6) 

where  £ is  can  be  any  chosen  x- coordinate  at  time  t-0.  The  integration 

of  Eq.  (B.5)  gives 

*u 

du  - 0 

.u(£) 

U - u(C)  (B . 7) 

where  u(£)  is  the  solution  at  the  coordinate  x-£  and  time  t-0. 

Equation  (B.7)  shows  that  along  each  characteristic  curve  given  by 
Eq . (B . 6)  the  solution  u does  not  change  and  can  be  determined  once  the 
initial  and  boundary  conditions  of  the  problem  are  known. 

For  the  problem  described  in  Section  7.1,  the  following  initial  and 
boundary  conditions  were  given: 
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u(x, 0) 


2(4  - x) 
3 
1 


1 < x < 2.5 


x > 2.5 


(B.8) 


u(Lltt)  - 2 (B . 9) 

In  the  equations  above  x varies  between  L^-l  m and  L2“ll  m,  and  the 
constant  in  Eq . (b.l)  is  c - 1 m/sec. 

The  solution  for  this  problem  can  be  obtained  along  each  character- 
istic curve  given  by  Eq.  (B.6)  and  a chosen  x-coordinate  £.  Let  £-1, 
for  example.  Then,  the  characteristic  curve  associated  to  this  point  is 

x - 1 + t (B.10) 

and  the  solution  given  by  Eq.  (B.8)  at  time  t-0  and  x-1  is 

u(0,l)  - 2 (B.ll) 

According  to  Eq.  (B.7)  this  solution  is  valid  along  the  character- 
istic curve  given  by  Eq.  (B.10).  Thus,  for  t]_-1.5  sec,  t2-  3.0  sec, 
t3«4.5  sec,  t4~6.0  sec,  t5=7 . 5 sec  and  tg-SLO  sec,  Eq.  (B.10)  gives  the 
correspondent  coordinates  where  the  solution  is  the  same  u-2 . These 
coordinates  are  x^-2.5  m,  X2-4.0  m,  X3-5.5  m,  X4-7.O  m,  X5-8 . 5 m and 
X6=10  m.  The  same  procedure  is  repeated  to  determine  the  solution 
along  any  other  characteristic  curve.  The  results  are  shown  in  Fig.  7.1 
and  represent  the  exact  solution  for  the  linearized  inviscid  Burgers' 
equation. 


APPENDIX  C 


THE  SAAGG  METHOD  AND  COMPUTATIONAL  EFFICIENCY 


The  SAAGG  method  can  greatly  reduce  the  number  of  grid  points  needed 
to  obtain  solutions  of  partial  differential  equations  by  using  finite- 
difference  methods  (FDMs).  This  reduction  in  the  number  of  grid  points 
enhances  the  computational  efficiency  of  the  FDMs  when  the  SAAGG  method 
is  used  for  generating  non-equally  grid  points  in  the  physical  domain. 

Computational  efficiency  can  be  estimated  by  quantifying  the  amount 
of  computer  memory  used  and  the  processing  time  (i.e.  CPU  and  I/O  time) 
needed  to  complete  a computing  task.  In  this  appendix,  the  improvements 
in  computational  efficiency  gained  by  using  the  SAAGG  method  are  des- 
cribed. 


One- dimensional  Problems 

In  this  section,  the  processing  time  and  the  computer  memory  used  to 
obtain  two  solutions  of  the  one -dimensional  problem  presented  in  Section 
6.3  are  compared.  One  solution  was  obtained  by  using  a FDM  with  21  non- 
equally  spaced  grid  points  distributed  in  the  physical  domain  by  the 
SAAGG  method.  The  other  solution  was  obtained  by  using  the  same  FDM 
with  251  equally  spaced  grid  points  distributed  in  the  physical  domain 
so  that  the  minimum  grid  spacing  was  equal  to  the  minimum  grid  spacing 
obtained  by  using  the  SAAGG  method. 
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The  exact  solution  of  that  problem  and  the  solutions  obtained  by 
using  non-equally  spaced  and  equally  spaced  grid  points  are  shown  in 
Fig.  C.l.  This  figure  shows  that  in  order  to  obtain  solutions  with 
approximately  the  same  accuracy  it  was  necessary  to  use  251  equally 
spaced  grid  points  instead  of  the  21  non-equally  spaced  grid  points 
generated  by  using  the  SAAGG  method.  This  means  that  the  number  of  grid 
points  was  increased  by  a factor  of  12  and  consequently  the  computa- 
tional efficiency  was  decreased.  The  processing  time  and  the  computer 
memory  needed  to  obtain  these  two  solutions  are  compared  in  Table  C.l. 
From  this  table  it  can  be  seen  that  by  using  the  SAAGG  method  instead  of 
using  equally  spaced  grid  points  the  processing  time  and  the  computer 
memory  can  be  reduced  by  a factor  of  2.33  and  4.28,  respectively.  Here, 
it  is  noted  that  the  problem  described  in  Section  6.3  is  very  simple  and 
does  not  require  much  computational  efforts  to  be  solved.  Hence,  the 
ratio  between  the  increasing  in  the  processing  time  per  grid  point  for 
using  the  SAAGG  method  (0.56-0.11)  and  the  processing  time  per  grid 
point  for  solving  the  problem  without  the  SAAGG  method  (0.11)  is  ap- 
proximately equal  to  4.0.  This  ratio  is  reduced  by  orders  of  magnitude 
when  the  problem  to  be  solved  is  more  complex  because  in  this  case  the 
processing  time  expended  just  for  using  the  SAAGG  method  will  be  a small 
fraction  of  the  total  processing  time. 

Multidimensional  Problems 

For  multidimensional  problems,  two  cases  can  be  considered.  The 
first  is  when  the  solutions  are  obtained  by  using  explicit  FDMs  and  the 
second  is  when  the  solutions  are  obtained  by  using  implicit  FDM. 
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Figure  C.l  Solution  obtained  for  the  initial-value  problem  given 
by  Eqs . (6.3.1)  and  (6.3.2).  The  units  for  x and  u are 
m and  m/sec,  respectively. 

(A)  Numerical  solution  obtained  with  251  equally  spaced 
grid  points. 

(B)  Numerical  solution  obtained  with  21  non-equally 
spaced  grid  points  generated  by  using  the  SAAGG 
method. 

(C)  Exact  solution 
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Table  C.l* 


Comparison  of  Computational  Efficiency 


SAAGG  METHOD 

EQUALLY  SPACED  METHOD 

Number  of  Grid  Points 

21 

251 

Normalized  Minimum  Grid 
Points  Axmin/(L2-L1) 

0.004 

0.004 

Processing  Time  Tt  (sec) 

11.75 

27.46 

Processing  Time  Per  Grid 
Point  Tp  (sec/point) 

0.56 

0.11 

Memory  (Byte) 

105x8 

502x8 

* All  computing  done  on  a Zenith  PC-ZW-158-43 
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Explicit  Methods 

When  using  explicit  FDMs  to  obtain  solutions  for  multidimensional 
problems,  the  processing  time  is  approximately  equal  to  the  processing 
time  per  grid  obtained  for  one -dimensional  problems  multiplied  by  the 
number  of  dimensions  (i.e.  2 for  2-D  and  3 for  3-D)  and  by  the  total 
number  of  grid  points.  Thus,  the  processing  time  for  a multidimensional 
problem  can  be  expressed  by 

Tt  - D Tp  (ILxJL)  (C.l) 

where  Tt  and  Tp  are  defined  in  table  C.l,  D is  the  number  of  dimensions 
and  IL  and  JL  are  the  number  of  points  in  each  grid  line  in  the  x-  and 
y-direction,  respectively. 

From  the  above  equation  it  can  be  seen  that  if  the  number  of  grid 
points  in  x-  and  y-direction  were  multiplied  by  a and  b,  respectively, 
the  time  processing  would  be  multiplied  by  a factor  equal  to  a times  b. 

Taking  for  example  the  same  figures  shown  in  Table  C.l  the  proces- 
sing time  for  a two-dimensional  problem  would  be  8 min  and  14  sec  by 
using  the  SAAGG  method  and  3 h and  51  min  by  using  equally  spaced  grid 
points.  For  a three-dimensional  problem,  the  time  processing  would  be 
12  min  and  21  sec  by  using  the  SAAGG  method  and  5 h and  21  min  by  using 
equally  spaced  grid  points. 

Implicit  Methods 

In  order  to  solve  a two-dimensional  problem  by  using  implicit  FDMs, 
the  following  matrix  equation  needed  to  be  solved: 
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[A]x  - b (C . 2) 

where  the  matrix  [A]  is  of  order  N and  N-ILxJL. 

The  number  of  arithmetical  operations  necessary  for  solving  Eq.(C.2) 
by  using  Gauss  elimination  is  of  order  of  0(N3) . Thus  if  the  number  of 
grid  points  in  x-  and  y-direction  were  multiplied  by  a constant  a and  b, 
respectively,  the  number  of  arithmetical  operations  would  be  multiplied 
by  a factor  (ab)3.  In  this  case  the  SAAGG  method  is  tremendously 
advantageous  because  by  reducing  the  number  of  grid  points  in  each 
coordinate  direction  by  order  of  10,  for  example,  the  number  of 
arithmetical  operations  will  be  reduced  by  order  of  10^. 
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